实验三、随机信号的功率谱估计方法

随机信号的功率谱估计方法 一、 实验目的 1、 利用自相关函数法和周期图法实现对随机信号的功率谱估计 2、 观察数据长度、自相关序列长度、信噪比、窗函数、平均次 数等谱估计的 分辨率、稳定性、主瓣宽度和旁瓣效应的影响。 3、 学习使用 FFT 提高谱估计的运算速度。 4、 体会非参数化功率谱估计方法的优缺点。 二、 实验原理 假设信号 x(n)为平稳随机过程，其自相关函数定义为 (3-（ ） ()(+) 1) 其中 E表示取数学期望，*表示共轭运算。根据定义， x(n)的功率谱密度 与()P 自相关序列 存在下面关系：()m (3-2)()()jmmPe (3-3)12jd 但是，实际中我们很难得到准确的自相关序列 ，只能通过随机信号的() 一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。 目前，常用的线性谱估计方法有两种，即相关函数法和周期图方法，本实 验对这两种方法分别予以讨论。 1 自相关函数法 假设我们已知随机信号x(n)的M长的自相关序列 ，利用自相关函数法()m 可以得到x(n)的功率谱估计： (3-4)*1()() Lmixi1()()MjmmPe 利用窗函数，上式又可表达为 (3-5)()()RjMmPWe 其中， 为矩形窗函数，定义为()RMW (3-6)1()0RM 因此， 实际上是真正功率谱 与窗函数 傅立叶变换的卷积。()P()P()RMWm 矩形窗函数不仅降低了谱估计的分辨率，而且使谱估计产生了旁瓣。为了降低 旁瓣影响，可以采用具有较小旁瓣的窗函数，如Hamming 窗，它定义为 (3-7)0.546cos()HMWm 这种窗函数可以有效的抑制旁瓣，但是，此时主瓣宽度增大，从而降低了 谱估计的分辨率，这种主瓣和旁瓣之间的矛盾在线性谱估计方法中是无法解决 的。 2 周期图方法 假设已知随机信号 的N个样本，利用周期图方法，信号x(n)的功率谱估()xn 计为 (3-8) 210()()NjnnPxe 利用上述方法得到的谱估计方差与信号 的功率谱平方 成正比，()2()P 为了减小它的方差，可以将信号序列进行分段处理，然后再求各分段结果的平 均，这就是平均周期图方法，即Bartlett方法。 （1）Bartlett 平均周期图方法 将一个随机序列 (0nN）分成K段，每段长度为L，各段之间互不重()xn 迭，因而N=KL，可以想到，第i段的信号序列可表示为 (3-9)()(1)ii0,1.,niK 对于每一段的周期图又可写成 , (3-10) 210()()LjniinIxe1.,i 于是，功率谱估计定义为 (3-11)1()() KiPI 因此，对于固定的记录长度来讲，分段数K增大可使谱估计的方差减小，但 是由于L的减小，相应的功率谱主瓣增宽，谱分辨率降低，显然，方差和分辨率 也是矛盾的。 除了分辨率降低以外，分段处理还会引起序列的长度有限所带来的旁瓣效 应。为减小这种影响，最有效的办法是给分段序列用适当的窗函数加权，可以 得到较平滑的谱估计，当然，相应的分辨率也有所下降。 （2）平滑平均周期图方法 这时一种改进的Bartlett周期图方法，它特别适用于FFT直接计算功率谱估 值。将长度为N的平稳随机信号序列x(n)分成K段，每段长度为L，即L=N/K。但 这里在计算周期图之前，先用窗函数 给每段序列 加权，K个修正的周()LWn()ixn 期图定义为 , (3-12) 210()()Lwjni iLnIxeU1.,i 其中U表示窗函数序列 的能量，()LW (3-13)201() Ln 在这种情况下，功率谱估计可按下面表达式给出： (3-14)1()() KwiPI 本实验主要是利用自相关函数法和周期图方法对下面受噪声干扰扰的正弦 信号进行谱估计： (3-15)()1()(i NSjnixae 其中NS为正弦个数， 分别为第i个正弦信号的数字频率、相位和幅度，,和 随机的分布在（ 0, 2）之间，w(n)为零均值方差等于 的复高斯白噪声。 2 三、实验内容和步骤 1 仔细阅读有关线性谱估计的内容，根据给出的框图3-1编制自相关函数法谱 估计的程序。运行程序，输入N=100,M=10, 选择 1=0.6,1=0,1=1, 2=0 矩形窗。观察谱峰位置是否正确（注意：由于窗效应可能引起谱估计的非正定） 。 如下图所示， 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 14 16 18 20 估估估估估估估估 图（1） 由上图（1）的仿真结果看出，谱峰的位置处于0.6 处，故估计较精确。 开始 输入参数：数据长度 N， 自相关函数个数 M，平均 次数 K 信号产生：输入正弦个数 Ns，每个正弦信号的数字频率、 相位和幅度，白噪声信号的方差 ，按照公式（3.16） 2 产生复正弦加白噪声信号的 N 个采样 结束 自相关函数法 由 N 个 x(n)估计出自相关 序列(M 长),并对此自相关 序列加矩形窗或 Hamming 窗，利用公式(3.5)计算 0，2）之间的 128 个功 率谱抽样点 周期图方法 输入 FFT 点数 ，按照 公式(3.11 )、(3.12 ) 和(3.13)计算 NF点功率 谱 图 3-1 2 观察并记录参数变化对谱估计性能的影响。 （1）改变M=5，其它输入同步骤1，观察功率谱估计的主瓣宽度和旁瓣大小随自 相关序列长度的变化情况。 0 0.2 0.4 0.6 0.8 10 1 2 3 4 5 6 7 8 9 估估估估估估估估 图（2） 上图（2）为 M=5 是，矩形窗函数的功率谱估计的仿真图。比较 M=5 与 M=10 时的图形可以看出：M 越小，主瓣宽度越大，分辨率越低，谱峰高度越低 低。 （2）选择窗函数为 Hamming 窗，其它输入同步骤 1，观察不同的窗函数对谱估 计性能的影响。 0 0.2 0.4 0.6 0.8 10 1 2 3 4 5 6 7 8 9 10 hamming估估估估估估 图（3） 比较图（1）图（3）可以看出，图（3）的主瓣宽度增加了近一倍，旁瓣相 应也减少了许多，同时发现 Hamming 窗的分辨率降低了。 （3）改变 ，其它输入同步骤 1，观察初始相位的变化对谱估计性能的4/1 影响。 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 14 16 18 20 估估估估估估估估 图（4） 比较图（1）、（4）可以看出，说明相位变化对谱估计的性能没有影响。 （4）改变 ，其它输入同步骤 1，观察信噪比变化对谱估计性能的影响。 2=1 （5）改变 N=10， ，其它输入同步骤 1，结合(4)的内容，观察数据长度2w 及信噪比对谱估计性能的影响。 0 0.2 0.4 0.6 0.8 10 5 10 15 20 25 估估估估估估估估 0 0.2 0.4 0.6 0.8 10 5 10 15 20 25 估估估估估估估估 图（5） 图（6） 如上图（5）为 时，加矩形窗普估计的仿真图，对比图（1）可以看 2=1 出，当增大高斯白噪声的方差时，对谱估计的峰值没有产生影响，但对旁瓣的 宽度和幅值有影响。 如上图（6）所示，在 N 较大时，改变信噪比对谱估计影响比较小，而当 N 较小时，改变信噪比对谱估计的影响较大。有图（6）可以看出此时频谱峰值已 偏离了 0.6 ，同时旁瓣的幅值增大了许多。 3运行程序，输入 ，1,0,8.0,6.,2,10 212121 aS ，选择矩形窗，调整自相关序列长度M，使得两个正弦频率分量临界分辨02w 出来，纪录此时的M值，并绘制此时的功率谱图。同样，在加Hamming窗的情况 下，记录使两个正弦频率分量临界分开的M值，并绘制此时的功率谱图。 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 估估估估估估估估 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 估估估估估估估估 图（7） 图（8） 图（7）、（8）分别是M=7和M=8时，加矩形窗的普估计仿真图，可以看出 M=8时可以使两个正弦频率分量临界分辨出来。 0 0.2 0.4 0.6 0.8 10 12 34 56 78 910 haming估估估估估估 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 haming估估估估估估 图（9） 图（10） 图（9）、（10）分别是M=9和M=10时，加Hamming窗的普估计仿真图，可以看出 M=10时可以使两个正弦频率分量临界分辨出来。 4运行自相关函数法谱估计程序，输入 ，选择矩210,0,1wNMS 形窗，观察利用自相关函数法得到的白噪声信号谱估计。改变 M=3，20，观察 M 的变化对白噪声谱估计的影响。 如下图（11）、（12）、（13）分别为 M=3、10、20 时，加矩形窗的普估 计的仿真图，可以看出，M 越大，谱估计的各个分量越越清晰。 0 0.2 0.4 0.6 0.8 10 1 2 3 4 5 6 估估估估估估估估 图（11） 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 14 16 18 估估估估估估估估 0 0.2 0.4 0.6 0.8 10 5 10 15 20 25 30 35 估估估估估估估估 图（12） 图（13） 5. 根据框图，编制周期图法谱估计程序。自程序 FFT 可以直接调用 Matlab 中 的函数。运行程序，输入 。选择窗函数为矩形窗,NFK ， ，观察谱峰位置是否正确，并与步骤 1 结1NS210.6,0wwa 果比较。 0 0.2 0.4 0.6 0.8 10 1 23 45 6 78 910 估估估估估估估 图（14） 对比图（1），上图中的谱峰位置同样也出现在0.6 处，谱的幅值明显降低， 主瓣宽度增加，旁瓣数量减少。 6利用周期图方法重复步骤2、3、4的内容，这里L=N/K相当于自相关函数中的 M，观察周期图法谱估计和自相关函数法谱估计在分辨率和稳定性方面的差别。 1）、 0 0.2 0.4 0.6 0.8 10 0.51 1.52 2.53 3.54 4.55 估估估估估估估 0 0.2 0.4 0.6 0.8 10 1 2 3 4 5 6 7 8x 10-3 估haming估估估估 图（15） 图（16） 如图（15）、（16）分别为N=100，K=20时，加矩形窗和汉明窗后的普估计 图。N/K减小，周期图法谱估计的主瓣宽度越大，分辨率降低，幅值降低，而且 旁瓣数量减少。这和自相关函数法中减小M时，谱估计出现的规律是一致的。 加汉明窗时的普估计，对比于图（16）和图（3）可以看出，周期图法得到 的谱的主瓣宽度加倍，旁瓣减小许多。周期图法谱估计和自相关法谱估计出现 的规律相似。 2）、改变 ，其它输入同步骤 1，观察初始相位的变化对谱估计性能的4/1 影响。 0 0.2 0.4 0.6 0.8 10 1 2 3 4 5 6 7 8 9 10 估估估估估估估 图（17） 对比图（14）以看出，改变初始相位的值对谱估计没影响。 改变 ，其它输入同步骤 1，观察信噪比变化对谱估计性能的影响。 2=1 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 估估估估估估估 图（18） 如上图所示，增大高斯白噪声的方差时，对谱估计的峰值没有影响，但却 使旁瓣增大。 3）、改变 N=50， ，其它输入同步骤 1，结合(4)的内容，观察数据长度12w 及信噪比对谱估计性能的影响。 0 0.2 0.4 0.6 0.8 10 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 估估估估估估估 图（19） 由图可知，当N减小，同时增大噪声方差时，矩形窗的周期图谱估计的主瓣 宽度增加。 4）、运行程序，输入 ，1,0,8.0,6.,2,10 212121 aNS ，选择矩形窗，调整自相关序列长度M，使得两个正弦频率分量临界分辨02w 出来，纪录此时的M值，并绘制此时的功率谱图。同样，在加Hamming窗的情况 下，记录使两个正弦频率分量临界分开的M值，并绘制此时的功率谱图。 固定K=10，分别取N=30和40，运行程序可以得到下图（20）、（21）所示 的加矩形窗普估计的图形。 0 0.2 0.4 0.6 0.8 10 1 2 3 4 5 6x 10-15 估估估估估估估 0 0.2 0.4 0.6 0.8 10 0.5 1 1.5 2 2.5x 10-15 估估估估估估估 图（20） 图（21） 图（20）、（21）分别是N/K=3、4时，加矩形窗的普估计仿真图，可以看 出N/K=4时可以使两个正弦频率分量临界分辨出来。只不过此时谱峰的位置已不 在0.6 和0.8 处了，普估计不准确了。 5）、固定K=10，分别取N=80和90，运行程序可以得到下图（22）、（23）所示 的加汉明窗普估计的图形。 0 0.2 0.4 0.6 0.8 10 0.5 1 1.5 2 2.5 x 10-17 估haming估估估估 0 0.2 0.4 0.6 0.8 10 1 2 3 4 5 6x 10-17 估haming估估估估 图（22） 图（23） 图（22）、（23）分别是N/K=8、9时，加汉明窗时的普估计仿真图，可以 看出N/K=9时可以使两个正弦频率分量临界分辨出来。此时估计仍发生了错误。 6）、运行自相关函数法谱估计程序，输入 ，选择210,0,1wNMS 矩形窗，观察利用自相关函数法得到的白噪声信号谱估计。改变M=3，20，观察 M的变化对白噪声谱估计的影响。 固定N=100，分别取K=50，10，5，即N/K=2、10、20，运行程序，对应的有 图（24）、（25）、（26）的仿真结果。 0 0.2 0.4 0.6 0.8 10 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 估估估估估估估 图（24） 0 0.2 0.4 0.6 0.8 10 2 4 6 8 10 12 14 估估估估估估估 0 0.2 0.4 0.6 0.8 10 5 10 15 20 25 估估估估估估估 图（25） 图（26） 对比以上三图可以看出，N/K的值越大，谱估计的各个分量就越清晰，也就 越容易分辨了。 四、实验主要结论 1、使用不同的窗函数谱估计的质量是不一样的，矩形窗的主瓣较窄，分辨率较 好，但方差较大，噪声水平较高；而 Hamming 窗的主瓣较宽，分辨率较低，但 方差较小，噪声水平较低。 2、周期图法在分辨率和稳定性方面是优于自相关函数法的谱估计：增大 M(或 N/K)值，会使主瓣宽度减小，即分辨率增加；反之，M 越小，主瓣越宽，分辨 率越小。另一方面，对于周期图法，仍存在着频率分辨率低、方差性能不好的 问题，原因是谱估计时对数据加窗截断，用有限个数据或其自相关函数来估计 无限个数据的功率谱，当数据很短时，这个问题更为突出。 3、初始相位对谱估计没有影响。 五、思考题 1. 证明：式（3-4）可以表示为： 10()2Re()(0)MjwmPwale 其中，Real表示取实部。 证明： ()= 1=+1() =1=0()+ 0=+1()(0) =1=0()+1=0()(0) 又自相关函数为偶函数，则上式可化简为： () =1=0()+1=0()(0)=1=0()（ + ） (0) 即： -1-=0()2Real(m)e()MjP 2、 已知实信号 的 N 个样本 ，可以定义 的估计如(n)x(),.,(),xxN()AP 下： ()= 1=+1() 其中 |*0()()Nmixi 试证明： 210()()NjnAnPxe 证明： ()= 1=+1() = 1=+11|1=0()(+) 令 m+n=k,上式为： =11=0() 1=0()() =11=0() 1=0() 即有： 210()()NjnAnPxe 附录一 clear clc N=100; M=10; Ns=1; w1=0.6*pi; pha1=0; a1=1; sigma_w2=0; L=256; wn=sqrt(sigma_w2)*(randn(N,1); xn=zeros(N,1); for n=1:N xn(n)=a1*exp(j*(w1*n+pha1)+wn(n); end r_xx=xcorr(xn,unbiased); win_p=2*M-1; W_rect=rectwin(win_p); rect1=r_xx(N-M+1:N+M-1).*W_rect; rect2=fft(rect1,L); rect3=abs(rect2(1:L/2); figure(1),clf stem(0:L/2-1)/(L/2),rect3,b)%归一化 title(矩形窗功率谱估计); W_han=hamming(win_p); han1=r_xx(N-M+1:N+M-1).* W_han; han2=fft(han1,L); han3=abs(han2(1:L/2); figure(2),clf stem(0:L/2-1)/(L/2),han3,b) title(hamming 窗功率谱估计); 附录二 clear clc N=100; M=10; Ns=2; w1=0.6*pi; w2=0.8*pi; pha1=0; pha2=0; a1=1; a2=1; sigma_w2=0; L=256; wn=sqrt(sigma_w2)*(randn(N,1); xn=zeros(N,1); for n=1:N xn(n)=a1*exp(j*(w1*n+pha1)+a2*exp(j*(w2*n+pha2)+wn(n); end r_xx=xcorr(xn,unbiased); win_p=2*M-1; W_rect=rectwin(win_p); rect1=r_xx(N-M+1:N+M-1).*W_rect; rect2=fft(rect1,L); rect3=abs(rect2(1:L/2); figure(1),clf stem(0:L/2-1)/(L/2),rect3