正弦信号频率误差估计

一、对00()cos(2)()s t a f t m t πθ=++进行频率误差估计 1.插值FFT 估计频率原理
单一频率实正弦信号表示为
)2cos()(00θπ+=t f a t s (1)
其中0f a 、 和0θ分别为正弦信号的幅度、频率和初相。

按等间隔N T t /=∆对)(t s 在0～T 区间内进行采样得到长度为N 的序列)(n s 。

)(n s 的N 点DFT 记为)(k S ，鉴于实序列的DFT 的对称性，忽略DFT 频谱的负频率成分，只考虑离散频谱的前 N/2点，有
12/,...,2,1,0]},)(1
[exp{]/)(sin[2)](sin[)(0000-=---∙--=
N k T f k N
N j N T f k T f k a k S πθππ(2)
)(k S 幅度最大值处的离散频率索引值记做1k ，]int[01T f k =，]int[x 表示取最接近
x 的整数，对于较大的N ，在幅度最大处，)(k S 的幅度可以近似表示为
πδ
πδ2)
sin()(11Na k S A =
= (3)
其中f f k f ∆∆-=/)(10δ为信号频率与其DFT 幅度最大处对应频率的相对偏差，
T f /1=∆，δ的变化范围为-0.5～0.5。

在紧邻1k 的左侧和右侧的两条谱线中幅度
较大处（以下称为幅度次大值,对应的离散频率索引值记做1,122±=k k k ），)(k S 的幅度可近似表示为
)
1(2)sin()(22δππδ-=
=Na k S A (4)
2A 与1A 的比值记做α,根据式(3)和(4)式,有
δ
δα-=
=112
A A (5) 根据2A 与1A 的比值可以得到δ的估计值
2
121A A A +=
+=
αα
δ (6) 根据δ值可对由离散频谱得到的0f 的估计值进行插值从而得到更精细的频率估计值
f k f ∆±=∧
)(10δ (7)
式中符号根据2k 的位置确定，若112+=k k 取加号，反之取减号。

2．高斯白嗓声背景中插值FFT 方法的频率估计误差
在加性噪声背景中的观侧值采样序列可以表示为)()()(n z n s n r +=，)(n z 为高斯白噪声序列，其方差为2z σ，采样序列的信噪比为)2/(220z a SNR σ=。

)(n r 的DFT 可表示为
12/,...,2,1,0),()()(-=+=N k k Z k S k R (8)
)(k Z 的方差为
2)](var[z N k Z σ= (9)
在)(0SNR N 较大的情况下,采样序列)(n r 的DFT 在最大谱线处的幅度可近似表示为
11111)(U A Z real A R +=+= (10)
1Z 为迭加在)(k S 幅度最大处的噪声, )(11Z real U =为1Z 的实部。

同理，当δ不接近0（即次大谱线的幅度不接近0）时，DFT 次大谱线的幅度可表示为
22222)(U A Z real A R +=+= (11)
2Z 为迭加在)(k S 幅度次大处的噪声, )(22Z real U =为2Z 的实部。

1Z 和2Z 均与)(k Z ，因此1U 和2U 均服从高斯分布，而且不加窗时1U 和2U 相互独立，
2/)var()var(221z N U U σ==。

因此噪声背景中采样序列的DFT 幅度次大值与最大值之比α可表示为
1
12
2121
2/1/1A U A U A A R R ++∙
=
=
α (12) 在信噪比较高或N 较大情况下，22/A U 及11/A U 大于或接近1的概率很小，在下面分析估计误差时忽略这种小概率情况而假设1/22〈〈A U 及1/11〈〈A U ，因此可以将上式展开成级数，并路去高次项，α可近似表示为
)1(1
12212A U
A U A A -+=
α (13)
因此噪声背景中δ的估计值可近似表示为（以下在不致引起混淆的情况下忽略δ的绝对值符号）
)
(111
111
22211
12
222
1U A A U A A U A A U A A A -++-+∙+=+=∧
ααδα
(14) 同理也可以将式(14)展开成级数，并略去高次项，于是的δ估计值可近似示为
2
11
222
2111222211212)()(1)()(U A A U A A U A A U A A A A A A -+--+++=
∧
αδ (15)
由于上式的最后一项远小于前两项，可以忽略不计，所以上式又可以近似为
)()(11
2
22
211212U A A U A A A A A A -+++=
∧
αδ (16) 所以，αδ∧
的方差可表示为
)
(4)(2)1()()var(02
42122212212242121SNR Na A A A A N A A A A A z ++=
++=∧
σδα (17) 将1A 和2A 代人式(17) ，整理，并考虑到∧∧=δδα)(E ，得αδ∧
的均方误差
)
(sin )(])1[()1()var()(2
0222δδδδδδααc SNR N mse +--=
=∧
∧ (18)
当δ=0.5时，αδ∧
的均方误差最小，为
)
(32)(02
SNR N mse πδα=
∧
(19)
当δ接近0时，上式分析不成立。

当δ=0时，02=A ，2/1Na A =，22Z R =，
11A R ≈，)/(2/212Na Z A Z =≈α。

当δ接近0时，δα≈，αδ∧
的均方误差可表示为
)(2
)
()()()()()(02
22
2
22
SNR N Na N Na Z E E mse mse z ===
==∧
σααδα (20) 而式(18)在δ=0处的值为)(/10SNR N 。

DFT 系数幅度最大值1R 及次大值2R 都位于sinc 函数的主瓣内，112±=k k 。

当δ接近0时，2R 较小，与最大谱线另一侧第一旁瓣内的谱线幅度3R 接近，3R 的位置记做3k ，223±=k k 。

可以证明，无噪声干扰时始终有32R R >。

但是，在噪声干扰背景下，当δ较小时会发生32R R <的情况。

根据前面的讨论，δ的符号根据第二大谱线是位于最大谱线的左边还是右边来确定，因此DFT 幅度第二大谱线位置的错误将造成δ的符号错误。

当信号的实际频率为f k f ∆+=)(10δ时，根据DFT 系数幅度次大值与最大值的比值得到的频率估计值为
f k f ∆-=∧
∧
)(10δ，从而造成频率估计出现()f ∆δ2的误差。

若这种错误出现的概率
较大，则这种频率估计方法失去意义。

发生这种情况的概率与DFT 点数N 以及信噪比有关，同时也与δ值本身的大小（即信号频率偏离最大谱线的大小）有关。

因此对应一定的DFT 点数和信噪比，存在一个阈值TH δ，只有在TH δδ>时，才能利用DFT 系数的幅度信息进行频率插值。

为了确定上述情况造成的δ估值的误差，需首先求出发生这种判断错误的概率。

根据式(2)，()S k 在3k 处的幅度为
3sin()2(1)
Na A πδπδ=
+ (21)
因此3R 的幅度可以表示为
33333()R A real Z A U =+=+ (22)
将3R 当作第二大值的概率为32()()e P P R R δ=>。

代入3R 及2R ，整理得
32())e P P U U δ=->
(23)
因23232()()0,var()var()/2z E U E U U U N σ====，又3U 与2U 均服从正态分布且互不相关，有
32)]1U U -= (24) 所以
1
()2e P erfc δ= (25)
其中,()erfc x 为补误差函数。

因此，给定DFT 长度Ｎ、信噪比和δ，利用上
式可计算出相应的()e P δ。

根据()e P δ可计算次大值方向错误造成的频率估计的（相对）均方误差为
22()(2)()2e e mse P erfc δδδδ∧
== (26)
因此插值FFT 方法估计频率的总的(相对)均方误差为
2()()e T
mse mse ασδδ∧
∧
=+ (27)
将(18)和(26)代入(27)，得
222202
2
(1)[(1)]()sin 2()T
erfc N SN R c δδδδσδ+--+= (28)
以下是在高斯白嗓声背景中用插值FFT 方法对00()cos(2)()
s t a f t m t πθ=++进行频率估计误差
图1 不加窗时频率估计均方根误差理论计算
SNR0从上到下为42、24、12dB
3. 加窗对频率估计误差的影响
为了抑制频谱泄漏，通常在进行FFT 之前对采样序列进行加窗处理。

加窗后1A 和2A 仍随δ变化，但变化规律不同于式(3)和式(4)，因此不能用式(6)估计δ值。

对于Hanning 窗，加窗后仍可得到利用1A 和2A 的比值α估计的δ解析式，而对于其它窗函数可以采用数值方法获得利用α估计δ的对应关系。

下面以Hanning 窗为例，讨论加窗情况下插值FFT 方法的频率估计误差，在加Hanning 窗情况下，可按下式估计δ
211
αδα-=+ (29)
将(13)式代入(29)式，得
2112
212
12121
23()()A A A A U U A A A A A δ-=
+-++ (30)
10
10
10
10
10
由于加窗后2U 和1U 相关， 不能像式(17)那样直接计算δ的方差，根据公式
22var()2x y xy x y x y σσρσσ±=+±，其中2
x
σ和2y σ分别为x 和y 的方差，xy ρ为x 和y 的相关系数，有
22
122
21
1212421211
9var()[var()var()2()()]()A A A U U std U std U A A A A δρ∧
=--+ (31) std(x)表示x 的标准差，窗函数记为w(n)，则有
1
2
2
2120
var()var()()N z
z n U U n G N ωσ
ω
σ-====∑ (32)
式中，120
()/N n G n N ωω-==∑。

将式(32)代入式(31)得
22
1212124
120
9(2)var()4()A A A A NG A A SNR ω
ρδ∧
+-=+ (33) 式中12ρ为1U 和2U 的相关系数，12ρ与窗函数有关，且可由下式确定
1
2
2(/)
1
2
2(/)0
12121
2
120
()cov(,)
1()()()
()
N j n N N j n N n N n n n e U U n e std U std U NG n ππω
ω
ρω
ω
--=-====
=
∑∑∑ (34)
由于加窗使得正弦信号的FFT 主瓣变宽，3A 也在主瓣之内，其幅度及随δ的变化率与不加窗时相比均大大增加。

又由于3U 与2U 相关，(经过计算可得，对于 Hanning 窗，230.17ρ=)，出现32R R >的概率大大降低，次大值方向错误造成的频率估计的误差可以忽略，因此加窗情况下，δ的估计误差由式(33)确定。

10
10
10
图2加窗时频率估计均方根误差理论计算SNR0从上到下为42、24、12dB
以上是在高斯白嗓声背景中加Hanning 窗方法对
00()cos(2)()s t a f t m t πθ=++进行频率估计误差
二、对00()cos(2())s t a f t m t πθ=++进行频率误差估计
考虑收到的信号由单一频率复正弦信号和复高斯白噪声组成，表示为
0(),0,1,2,...,1t i t u t x Ae t N ωθ++==- (35)
其中0A ω、 和θ分别为正弦信号的幅度、频率和初相。

t x 是由0()i t t Ae z ωθ++近似而来，t z 是均值零、方差为2z σ的复高斯白噪声，12t t t z z jz =+，1t z 和2t z 是均值零、方差为2/2z σ的实高斯白噪声，且互不相关。

则可知式(35)中的t u 是方差为22/2z A σ的零均值实高斯白噪声。

则相位有
0,0,1,2,..., 1.t t x t u t N ωθ∠=++=- (36)
只希望进行频率估计，并避免差分相位数据的相位展开
1,0,1,2,...,2t t t x x t N +∆=∠-∠=- (37)
将式(36)代入(37)得
01t t t u u ω+∆=+- (38)
0ω的最大可能性估量等价于对式（38）进行最小方差无偏估计得
100(1)(1)T J C ωω-=-- (39)
其中[]012T
N -= ，[]1111T
= ，C 是t 的协方差矩阵，则
101111
T T C C ω-∧
-=
（40）
对应的0ω∧
的均方误差为
01
1
()11
T Var C ω∧
-=
（41） 因是0ω∧
无偏估计，有[]011,,0,1,2,...,2t t t u u u u t N ω+=+=-=- 。

而为了估计
1C -，首先可知t 是一个带有方差为22/2z A σ的噪声的实移动平均过程，系数
011,1b b ==-，协方差函数的系数为
222201
2
2
2
2012
2
(0)()2(1)(1)22()0, 2.
z z z z c b b A A c c b b A
A
c k k σσσσ=
+=
=-=
=-
=≥
则协方差矩阵可表示为
22
210001210020002z C A σ-⎛⎫ ⎪-- ⎪
= ⎪ ⎪⎝⎭
也可表示为
[]22min(,),1,12z ij ij C i j i j N A N σ⎡
⎤=
-≤≤-⎢⎥⎣⎦
（42） 经过以上计算之后，有
22
1
2
222120
(1)116(1)116T
Z
N T t t
t Z
N N A C N N A
C σωσ---=-=-=
∑
其中
2
3
22(1)2112t N t N N N ω⎧⎫⎡⎤--⎪⎪⎢⎥⎪⎪=-⎨⎬⎢⎥-⎪⎪
⎢⎥⎣⎦⎪⎪⎩⎭
因此式（40）可表示为
2
00
N t t t ωω-∧
==∑
由于2
1N t t ω-==∑，而0ω∧
在高信噪比下是无偏估计，则频率估计的式子可写成
11t t t t t x x x x *++∆=∠-∠=∠∠ （43）
因此
2
010
N t t t t x x ωω-∧
*+==∠∠∑ （44）
对应的0ω∧
的均方误差由式（41）可改写得
02
22
6()(1)
z Var A N N ωσ∧
=
- （45）
由式(43)和式(44)可将0ω∧
重写为
11
022
00126
(1)(1)N N t t t t t x x N N N N ω--∧
===∠-∠--∑∑ (47)。