1、数字信号处理实验报告姓名： 学号： 日期：2015.12.141. 实验任务信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为 10dB，长度为 ，信N号频率分别为 和 ，初始相位 ，取 ， 取不同的数值：1f20212.01sfsf10.3， 0.25。 为采样率。s（1 ） 分别用 Levinson 递推法和 Burg 法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。（2 ） 当正弦信号相位、频率、信噪比改变后，上述谱估计的结果有何变化？并作分析说明。2. 原理分析2.1 现代谱估计中的参数建模根据参数模型来描述随机信号的方法，我们可以知道，如果能确定信号 的信号xn
2、模型，根据信号观测数据求出模型参数，系统函数用 表示，模型输入白噪声，其方zH差为 ，信号的功率谱用下式求出：2w22iwjwxeeP按照这种求功率谱的思路，功率谱估计可分为三个步骤：（1）选择合适的信号模型；（2）根据 有限的观测数据，或者它的有限个自相关函数的估计值，估计模型的nx参数；（3）计算墨香的输出功率谱。其中以（1） 、 （2）两步最为关键。按照模型的不同，谱估计的方法有许多种，它们共同的特点是对信号观测区以外的数据不假设为 0，而先根据信号观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。下面分析 AR 谱估计的两种方法：自相关法
3、列文森（Levenson）递推法和伯格（Burg）递推法。这两种方法均为已知信号观测数据，估计功率谱，两者共同特点是由信号观测数据求模型系数时采用信号预测误差最小的原则。对于长记录数据，这些方法的估计质量是相似的，但对于短记录数据，不同方法之间存在差别。2.2 自相关法列文森（Levenson）递推法自相关法的出发点是选择 AR 模型参数使预测误差功率最小，预测误差功率为 21211 npiin nxaxNe假设信号 的数据区在 范围，有 个预测系数， 个数据经过冲激响应xn0PN为 的滤波器，输出预测误差 的长度为 ，因此应用下式计算：0,1pia ep2101210 PNnipiPNn n
4、xax的长度长于数据的长度，上式中数据 的两端需补充零点，相当于对无穷长的信en号加窗处理，得到长度为 N 的数据。上式对系数 的实部和虚部求微分使预测误差功率最pia小，得到（1） 120120pxxx xpxxx xrrrrarprrr 式中自相关函数采用有偏自相关估计，即 1,2,1,01*0 pmrnxNmxmnx对比上式，可知式（1）即为已推导出的 Yule-Walker 方程，因此自相关法也是基于解Yule-Walker 方程的一种方法。但是直接解该方程，需要计算逆矩阵，不方便，因此，基于 Yule-Walker 方程中自相关矩阵的性质，导出 Levinson-Durbin 递推法
5、，这是一种高效的解方程的方法。Levinson-Durbin 算法首先由一阶 AR 模型开始：一阶 AR 模型 的 Yule-Walker 方程为1p211010xxra由该方程解出 10xra2211,xr然后令 ，以此类推，可以得到一般递推公式如下：2,34p 称为反射系数， 。 ，随着阶数增加，预测误差功率将减pk1pk22p少或不变。由 k=1 开始递推，递推到 k=p，依次得到各阶模型参数，2221112,ppaa AR 模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算 222 1/pjwjwjwix iipeHeae这种方法计算简单，但需要预先估计出信号自相关函数，实
6、际中只能按照信号的有限个观测数据估计自相关函数。当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂（在信号谱峰附近产生虚假谱线） ；如数据很长，估计自相关函数较准确，但计算量大，应适当选择数据长度。2.3 伯格（Burg）递推法Levinson-Durbin 递推法需要由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算 AR 模型参数。伯格递推法利用 Levinson-Durbin 递推公式，导出前向预测误差与后向预测误差，并按照使它们最小的原则求出 ，从而实现不用估计自相关函数，直接用观测数据得出结果。pkBurg 递推法思想：借助格型预测误差滤波器，
7、求前向、后向预测误差平均功率，选择使其最小，求出 。之后，再利用 Levinson-Durbin 递推法求模型参数和输入噪声方pkpk差。设信号 的观测数据区间： ,前向、后向预测误差功率分别用 和xn01nNpf表示，预测误差平均功率用 表示，公式分别为pbp21Nfpfpne210Nbpbpne12pfpb前向、后向观测误差公式分别为 1pfpknxxnaepkbpxxe1*0pa上式中，信号项的自变量最大的是n,最小的是n-p，为了保证计算范围不超出给定的数据范围，在 和 计算公式中，选择求和范围为： 。pfb 1nN为求预测误差平均功率 最小时的反射系数 ，令 ,将前、后向预测误差的p
8、pk0p递推公式代入得 1212*NpnbpffnekBurg递推法求AR模型参数的递推公式总结如下：(1) 120,0Nx xnrxr (2) 0fbex,12N，(3) 1212*Npnbpffnek(4) 12ppk(5) ,iiiaa1,2ip，(6) ,pk(7) 1,2,111* Npneknefpbpbff 3. 编程思想（1 ） 编写程序产生题目要求的信号和噪声（2 ） 然后分别用两种方法的递推流程进行谱估计（3 ） 改变题目中要求的变量参数，分析结果的变化4. 代码Levensonclc;clear all;fs=100; %采样频率Ts=1/fs;N=27; %数据长度p1
9、=20; %阶数f1=0.2*fs;f2=0.25*fs; %设置信号频率pha1=0;pha2=0;%初始相位SNR=2; %设置信噪比%产生信号w=randn(1,N);Am=sqrt(2*10(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);xx=x1+x2;x=w+x1+x2;figure(1)subplot(2,1,1);plot(xx);ylim(-20,20);title(两个正弦信号波形 );grid on;subplot(2,1,2);plot(x);ylim(-20,20);ti
10、tle(加噪信号波形);grid on;%计算自相关函数R=;for m=1:Ns=0;for n=1:N-(m-1)s=s+x(n)*x(n+m-1);endR(m)=s/N;end %列文森递推法a(1,1)=-R(2)/R(1); cov(1)=R(1)+a(1,1)*R(2);for p=2:p1sum=0;for k=1:p-1sum=sum+a(p-1,k)*R(p-k+1);enda(p,p)=-(R(p+1)+sum)/cov(p-1);%计算反射系数for k=1:p-1a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k); endcov(p)=(1-(abs(
11、a(k,k)2)*cov(p-1);end %计算功率谱，功率谱以 2*pi 为周期，又信号为实信号，只需输出 0 到 P1 即可；W=0.01:0.01:pi;Z=0;for k=1:p1Z=Z+(a(p1,k).*exp(-j*k*W);endPSD=1./(abs(1+Z).2); %算出功率谱F=W*fs/(pi*2); %角频率坐标换算成频率坐标figure(2)plot(F,abs(PSD);xlabel(频率(Hz);ylabel( 功率谱);title(自相关法-列文森 Levenson 递推法的功率谱估计);grid;figure(3)p=1:p1;plot(p,cov(p)
12、,b);xlabel(模型阶数);ylabel(预测误差功率大小 );title(预测误差功率变化趋势 );grid;Burgclc;clear all;fs=100; %采样频率Ts=1/fs;N=28; %数据长度p1=20; %阶数 f1=0.2*fs;f2=0.25*fs; %设置信号频率pha1=0;pha2=0;SNR=2; %设置信噪比%产生信号w=randn(1,N); % 噪声为高斯白噪声，功率为 1.Am=sqrt(2*10(SNR/10);k=1:N;x1=Am*sin(2*pi*f1*k*Ts+pha1);x2=Am*sin(2*pi*f2*k*Ts+pha2);x=w
13、+x1+x2;%递推ef=zeros(p1,N);eb=zeros(p1,N);ef(1,:)=x; eb(1,:)=x; cov(1)=x*x/N; k(1)=0; a=zeros(p1+1,p1+1);for p=2:p1+1numerator=0;denominator=0;for n=p:Nnumerator= numerator+ (-2)*ef(p-1,n)*eb(p-1,n-1);denominator=denominator+(ef(p-1,n)2+(eb(p-1,n-1)2;endk(p)=numerator./(denominator+0.0001); cov(p)=(1-
14、abs(k(p)2).*cov(p-1);a(p,p)=k(p);for n=p:Nef(p,n)=ef(p-1,n)+k(p).*eb(p-1,n-1);eb(p,n)=eb(p-1,n-1)+k(p).*ef(p-1,n);endendk=k(1,2:p1+1);a=a(2:p1+1,2:p1+1);for p=2:p1for i=1:p-1a(p,i)=a(p-1,i)+k(p)*a(p-1,p-i);endend%功率谱Z=0;W=1:0.01:pi;for i=1:pZ=Z+(a(p,i)*exp(-j*W*i);end;pxx=cov(p1+1)./(abs(1+Z).2);F=W
15、*fs/(pi*2); %角频率坐标换算成频率坐标figure(1)plot(F,pxx);xlabel(Frequency(Hz);ylabel(Power Spectral Density);title(伯格（Burg）递推法的功率谱估计 );grid;% figure(2)% p=1:p1;% plot(p,cov(p),b);% xlabel(模型阶次 );% ylabel(预测误差功率大小);% title(预测误差功率的变化趋势);% grid;5. 实验结果及分析5.1 产生信号0 20 40 60 80 100 120 140-20-1001020 两两两两两两两两0 20 4
16、0 60 80 100 120 140-20-1001020 两两两两两两5.2 Levenson 递推法1. 数据长度的影响（阶数设置为 20 阶）N=320 5 10 15 20 25 30 35 40 45 50051015202530354045两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两N=640 5 10 15 20 25 30 35 40 45 50020406080100120140160180两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两N=1280 5 10 15 20 25 30 35 40 45 50010020030040
17、0500600两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两N=10240 5 10 15 20 25 30 35 40 45 50050010001500200025003000两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两N=81920 5 10 15 20 25 30 35 40 45 500100020003000400050006000两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两N=163840 5 10 15 20 25 30 35 40 45 500100020003000400050006000两两(Hz)两两
18、两两两两两-两两两Levenson两两两两两两两两两2. 模型阶数的影响（数据长度设置为 N=1024）预测误差功率的变化（从 1 阶到 50 阶）0 5 10 15 20 25 30 35 40 45 500510152025两两两两两两两两两两两两两两两两两两两两两两不同阶数（p1 表示）功率谱的图像如下：P1=50 5 10 15 20 25 30 35 40 45 50010203040506070两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两P1=100 5 10 15 20 25 30 35 40 45 50050100150200250300两两(Hz)两两
19、两两两两两-两两两Levenson两两两两两两两两两P1=200 5 10 15 20 25 30 35 40 45 50050010001500两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两P1=300 5 10 15 20 25 30 35 40 45 50050010001500200025003000350040004500两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两P1=500 5 10 15 20 25 30 35 40 45 500500100015002000250030003500两两(Hz)两两两两两两两-两两两Levenson
20、两两两两两两两两两3. 相位的影响（设置 N=128 p1=20）00 5 10 15 20 25 30 35 40 45 500100200300400500600两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两Pi/60 5 10 15 20 25 30 35 40 45 50050100150200250300350400450两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两Pi/30 5 10 15 20 25 30 35 40 45 50050100150200250300350400两两(Hz)两两两两两两两-两两两Levenson两两两两两
21、两两两两Pi/20 5 10 15 20 25 30 35 40 45 50050100150200250300350400450两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两Pi0 5 10 15 20 25 30 35 40 45 50050100150200250300350400450两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两4. 频率的影响（N=128 p1=20 初始相位为 0）Fs=100, f1=0.2*fs f2=0.25*fs0 5 10 15 20 25 30 35 40 45 50010020030040050060070
22、0两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两Fs=100, f1=0.2*fs f2=0.3*fs0 5 10 15 20 25 30 35 40 45 50050100150200250300350400两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两Fs=1000,f1=0.2*fs f2=0.25*fs0 50 100 150 200 250 300 350 400 450 5000100200300400500600两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两Fs=1000,f1=0.2*fs f2=0.3*fs0
23、50 100 150 200 250 300 350 400 450 500050100150200250300350两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两5. 信噪比的影响SNR =20dB 0 5 10 15 20 25 30 35 40 45 500200400600800100012001400两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两SNR =10dB0 5 10 15 20 25 30 35 40 45 500100200300400500600两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两SNR =6d
24、B0 5 10 15 20 25 30 35 40 45 50050100150200250300两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两SNR=2dB0 5 10 15 20 25 30 35 40 45 50020406080100120两两(Hz)两两两两两两两-两两两Levenson两两两两两两两两两5.3 Burg 递推法1. 数据长度的影响N=3215 20 25 30 35 40 45 5001002003004005006007008009001000Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两
25、两两N=6415 20 25 30 35 40 45 50050100150200250300350400450500Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两N=25615 20 25 30 35 40 45 500100200300400500600700Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两N=102415 20 25 30 35 40 45 50020004000600080001000012000Frequency(Hz)Power Spectral Den
26、sity两两两Burg两两两两两两两两两两N=819215 20 25 30 35 40 45 50050010001500200025003000350040004500Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两N=1638415 20 25 30 35 40 45 500100020003000400050006000Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两2. 模型阶数的影响（N=128）预测误差功率的变化趋势0 5 10 15 20 25 30 35 40 45
27、 500510152025两两两两两两两两两两两两两两两两两两两两两两两不同阶数的功率谱P1=515 20 25 30 35 40 45 50050100150200250Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两P1=1015 20 25 30 35 40 45 500100200300400500600700800Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两P1=2015 20 25 30 35 40 45 50050100150200250300350400Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两P1=5015 20 25 30 35 40 45 50050100150200250Frequency(Hz)Power Spectral Density两两两Burg两两两两两两两两两两
