一、 频偏估计
出现载波频率偏移（CFO），其原因通常有
发送端和接收端振荡器之间的频率不匹配多普勒效应
1.1 信号加频偏
若输入信号为
u
(
t
)
u(t)
u(t) ，则增加频偏后的信号
y
(
t
)
y(t)
y(t) 表示为
y
(
t
)
=
u
(
t
)
∗
(
c
o
s
(
2
π
∫
0
t
f
(
τ
)
d
τ
+
ψ
(
t
)
)
+
j
s
i
n
(
2
π
∫
0
t
f
(
τ
)
d
τ
+
ψ
(
t
)
)
)
y(t) = u(t)*\left(cos\left(2\pi \int_0^tf(\tau)d\tau + \psi(t)\right) + jsin\left(2\pi \int_0^tf(\tau)d\tau + \psi(t)\right)\right)
y(t)=u(t)∗(cos(2π∫0tf(τ)dτ+ψ(t))+jsin(2π∫0tf(τ)dτ+ψ(t))) 其中
f
(
t
)
f(t)
f(t) 表示频偏，
ψ
(
t
)
\psi(t)
ψ(t) 表示相偏
对于时域离散输入信号，输出可表示为
y
(
i
)
=
{
u
(
0
)
∗
(
c
o
s
(
ψ
(
0
)
)
+
j
s
i
n
(
ψ
(
0
)
)
)
i
=
0
;
u
(
i
)
∗
(
c
o
s
(
2
π
∑
n
=
0
i
−
1
f
(
n
)
Δ
t
+
ψ
(
i
)
)
+
j
s
i
n
(
2
π
∑
n
=
0
i
−
1
f
(
n
)
Δ
t
+
ψ
(
i
)
)
)
i
>
0
;
y(i) = \left\{ \begin{aligned} & u(0)*\left(cos(\psi(0))+jsin(\psi(0))\right) & i = 0; \\ & u(i)*\left(cos\left(2\pi \sum_{n=0}^{i-1}f(n)\Delta t + \psi(i)\right) + jsin\left(2\pi \sum_{n=0}^{i-1}f(n)\Delta t + \psi(i)\right)\right) & i > 0; \end{aligned} \right.
y(i)=⎩
⎨
⎧u(0)∗(cos(ψ(0))+jsin(ψ(0)))u(i)∗(cos(2πn=0∑i−1f(n)Δt+ψ(i))+jsin(2πn=0∑i−1f(n)Δt+ψ(i)))i=0;i>0;
其中
Δ
t
\Delta t
Δt 表示采样时长，其倒数
f
s
f_s
fs 为频域采样率，即
[
−
π
,
π
]
[-\pi, \pi]
[−π,π] 的频域分辨率
频偏最终都会累计体现在相位偏移上在星座图上，相偏会导致星座图旋转；对于每一个复数，相偏则改变其实部与虚部的大小对于某一时刻的信号，相偏为定值（假定系统只会产生相偏），故相偏会使信号星座图旋转一定的角度；而频偏对信号而言是一个累计量，更确切地说，频偏对系统是一个累积量，故同一个信号在不同时刻经过存在频偏的系统，输出信号的相位值是不一样的，而且进入系统的时间越晚，相位变化越大，体现在星座图上，也就是旋转的角度越大
1.2 频偏估计算法
基于相关性的频偏估计
文献 [1] 描述了一种基于相关性的估计算法，用于估计 PSK(Phase-shift keying) 与 PAM(Phase Amplitude Modulation) 信号的频偏。对
e
x
p
(
j
2
π
Δ
f
t
)
exp(j2\pi \Delta ft)
exp(j2πΔft) 求似然估计 (ML, Maximum Likelihood)，从而得到信号频偏
Δ
f
\Delta f
Δf
若观测到的信号为
r
(
k
)
=
e
j
(
2
π
Δ
f
k
T
s
+
θ
)
,
1
≤
k
≤
N
r(k)=e^{j(2\pi \Delta fkT_s+\theta)} , 1 ≤ k ≤ N
r(k)=ej(2πΔfkTs+θ),1≤k≤N 其中
T
s
T_s
Ts 为采样间隔，
θ
\theta
θ 为随机相位，
N
N
N 观测信号样本数量。利用似然估计计算频偏即寻找似然函数的最大值
Λ
(
Δ
f
)
≈
∣
∑
i
=
1
N
r
(
i
)
e
−
j
2
π
Δ
f
i
T
s
∣
2
=
∑
k
=
1
N
∑
m
=
1
N
r
(
k
)
r
(
m
)
∗
e
−
j
2
π
Δ
f
T
s
(
k
−
m
)
Λ(\Delta f)≈ \left|\sum^N_{i=1}r(i)e^{−j2\pi \Delta fiT_s}\right|^2=\sum^N_{k=1}\sum^N_{m=1}r(k)r(m)*e^{−j2\pi \Delta fT_s(k−m)}
Λ(Δf)≈
i=1∑Nr(i)e−j2πΔfiTs
2=k=1∑Nm=1∑Nr(k)r(m)∗e−j2πΔfTs(k−m) 化简后的目标函数可以看作是一个经过抛物线窗函数加权离散傅里叶变换，如下
I
m
{
∑
k
=
1
N
−
1
k
(
N
−
k
)
R
(
k
)
e
j
2
π
Δ
f
^
T
s
}
=
0
Im\left\{ \sum^{N−1}_{k=1}k(N−k)R(k)e^{j2\pi \Delta \widehat f T_s} \right\}=0
Im{k=1∑N−1k(N−k)R(k)ej2πΔf
Ts}=0
R
(
k
)
R(k)
R(k) 表示观测信号
r
(
k
)
r(k)
r(k) 估计自相关结果，定义为
R
(
k
)
:
=
1
N
−
k
∑
i
=
k
+
1
N
r
(
i
)
∗
r
(
i
−
k
)
,
0
≤
k
≤
N
−
1
R(k):=\frac1{N-k}\sum_{i=k+1}^Nr(i)*r(i-k) , 0 ≤ k ≤ N-1
R(k):=N−k1i=k+1∑Nr(i)∗r(i−k),0≤k≤N−1
k
(
N
−
k
)
k(N-k)
k(N−k) 项可以看作是一个抛物线窗函数。在文献 [1] 中也说明了当
k
=
0
k=0
k=0 或
k
k
k 接近
N
N
N 时，序列
r
(
k
)
r(k)
r(k) 的估计自相关结果
R
(
k
)
R(k)
R(k) 并不准确，同时在一定取值范围内，抛物线窗函数可以近似认为是值为 1 的矩形窗序列 ，进一步简化似然估计目标函数
I
m
{
∑
k
=
1
L
R
(
k
)
e
−
j
2
π
Δ
f
^
T
s
}
=
0
，
1
≤
L
≤
N
−
1
Im\left\{ \sum^L_{k=1}R(k)e^{-j2\pi \Delta \widehat f T_s} \right\}=0， 1 ≤ L ≤ N-1
Im{k=1∑LR(k)e−j2πΔf
Ts}=0，1≤L≤N−1 使上述目标函数成立，则频偏估计值
Δ
f
^
\Delta \widehat f
Δf
为
Δ
f
^
≌
f
s
a
m
p
π
(
L
+
1
)
a
r
g
{
∑
k
=
1
L
R
(
k
)
}
\Delta \widehat f≌\frac{f_{samp}}{\pi (L+1)}arg\left\{\sum_{k=1}^LR(k)\right\}
Δf
≌π(L+1)fsamparg{k=1∑LR(k)} 采样频率
f
s
a
m
p
f_{samp}
fsamp 是采样间隔
T
s
T_s
Ts 的倒数，且用于计算自相关序列的样本长度
L
L
L 定义为
L
=
r
o
u
n
d
(
f
s
a
m
p
f
m
a
x
)
−
1
L=round\left(\frac{f_{samp}}{f_{max}}\right)-1
L=round(fmaxfsamp)−1
f
m
a
x
f_{max}
fmax 是预测频偏的最大值， round 代表就近取整。当
L
≥
7
L ≥ 7
L≥7 且
f
m
a
x
≤
f
s
a
m
p
/
4
M
f_{max} ≤ f_{samp} / 4M
fmax≤fsamp/4M 时，该频偏估计算法可取得较好的性能
基于 fft 的频偏估计
基于 fft 的估计算法适用于估算不同调制方式下的信号频偏，一些粗频偏补偿器使用该算法以支持不同调制方式的信号处理
该算法在接收信号序列作
m
m
m 次乘方运算的基础上，使用周期图法估计频偏
Δ
f
^
\Delta \widehat f
Δf
，计算公式为
Δ
f
^
=
f
s
a
m
p
N
∗
m
a
r
g
m
a
x
f
∣
∑
k
=
0
N
−
1
r
m
(
k
)
e
−
j
2
π
k
t
/
N
∣
,
−
R
s
y
m
2
≤
f
≤
R
s
y
m
2
\Delta \widehat f=\frac{f_{samp}}{N*m}arg\;\underset{f}{max}\left|\sum_{k=0}^{N-1}r^m(k)e^{-j2\pi kt/N}\right|,-\frac{R_{sym}}2≤f≤\frac{R_{sym}}2
Δf
=N∗mfsampargfmax
k=0∑N−1rm(k)e−j2πkt/N
,−2Rsym≤f≤2Rsym 其中
m
m
m 为调制阶数，
r
(
k
)
r(k)
r(k) 为接收序列，
R
s
y
m
R_{sym}
Rsym 为符号速率，
N
N
N 为序列样本数量
实质上，就是在
[
−
R
s
y
m
/
2
,
R
s
y
m
/
2
]
[-R_{sym}/2, R_{sym}/2]
[−Rsym/2,Rsym/2] 范围内，找一个
f
f
f 使接收序列的
m
m
m 次方在时域的均值最大。该公式在形式上也可以看作是
r
m
(
t
)
r^m(t)
rm(t) 的离散傅里叶变换，寻找
f
f
f 使时域平均最大，等同于找
r
m
(
t
)
r^m(t)
rm(t) 频谱峰值所在的频点。
FFT 点数 N 计算方式如下
N
=
2
⌈
l
o
g
2
(
f
s
a
m
p
f
r
)
⌉
N=2^{\left\lceil log_2(\frac{f_{samp}}{f_r})\right\rceil}
N=2⌈log2(frfsamp)⌉ 其中
f
r
f_r
fr 为频率分辨率
文献 [1] Carrier frequency recovery in all-digital modems for burst-mode transmissions | IEEE Journals & Magazine | IEEE Xplore
通信中的同步（一）——基于FFT的信号信号载波频偏估计_8psk频偏估计-CSDN博客
Compensate for frequency offset of PAM, PSK, or QAM signal - MATLAB - MathWorks 中国