谱估计实验报告

1、第四章上机作业实验报告实验题目1、假设一平稳随机信号为，其中是均值为0，方差为1的白噪声，数据长度为1024。（1）、产生符合要求的和；（2）、给出信号x(n)的理想功率谱；（3）、编写周期图谱估计函数，估计数据长度N=1024及256时信号功率谱，分析估计效果。（4）、编写Bartlett平均周期图函数，估计当数据长度N=1024及256时，分段数L分别为2和8时信号的功率谱，分析估计效果。2、假设均值为0，方差为1的白噪声中混有两个正弦信号，该正弦信号的频率分别为100Hz和110Hz，信噪比分别为10dB和30dB，初始相位都为0，采样频率为1000Hz。（1）、采用自相关法、Burg法

2、、协方差法、修正协方差法估计功率谱，分析数据长度和模型阶次对估计结果的影响（可采用MATLAB自带的功率谱分析函数）。（2）、调整正弦信号信噪比，分析信噪比的降低对估计效果的影响。报告内容一、实验题目一1、问题分析（1）、w(n)与x(n)的产生w(n)产生：均值为0，方差为1白噪声利用matlab中randn函数即可。表达如下：w=sqrt(1)*randn(1,N); sqrt(1)表示方差为1。x(n)产生：第一种思路：利用迭代的方法由，其中，然后利用上述公式依次向后递推即可得。matlab代码实现如下，注意到matlab中元素下标都是从1开始的：x=w(1) zeros(1,N-1);

3、for i=2:N x(i)=0.8*x(i-1)+w(i);end此方法简单，可以很容易地产生所需数目的数据。第二种思路：利用卷积的方法对线性时不变系统，输入输出满足卷积关系：。由，可得，从而可得系统的冲击响应：。然后进行卷积运算即可。Matlab代码实现如下：n=1:N;h(n)=(0.8).(n-1); %要注意n-1不是n，因为冲击响应是从0开始的y=conv(w(n),h(n)； %共2*N-1项，取前N项即可需注意：实际h(n)是从0开始的，matlab处理元素从下标1开始，因此，公式中应是n-1不是n。而且，计算完成后卷积结果是为2*N-1项，取前N项即可。两种方法结果为方便观察

4、，令N=5时，实验结果如下：x = 0.6232 1.2976 1.9790 0.5911 0.6849 y = 0.6232 1.2976 1.9790 0.5911 0.6849 0.3437 0.0131 -0.2978 0.0868取卷积的前N项，可以看出两种方式结果是相同的。（2）、信号x(n)的理想功率谱系统为AR模型，理想功率谱为：因此，对h(n)进行傅里叶变换后，取模的平方即可。（3）、周期图法谱估计根据周期图法谱估计原理：先对观测数据x(n)进行傅里叶变换，后平方，最后除以N即可。（4）、Bartlett平均周期图谱估计为了减小估计的方差，将数据分为L段，则每段有M=N/L个

5、数据，分别用周期图法进行估计后求平均。具体公式如下：将得到的L个周期图进行平均，作为信号x(n)的功率谱估计，公式如下：经过平均处理，可使估计方差减小。2、实验结果与分析（1）、w(n)与x(n)的产生图1 白噪声w(n)图2 平稳随机信号x(n)（2）、信号x(n)的理想功率谱信号理想功率谱如下图3所示：图3 信号的理想功率谱从图中可以看出，理想功率谱是平滑的。下图4是功率谱的分贝形式：图4 信号的理想功率谱（dB）（3）、周期图法谱估计当数据长度N为1024时，实验结果如下图5所示：图5 N=1024周期图谱估计结果当数据长度N为256时，实验结果如下图6所示：图6 N=256周期图谱估计

6、结果对比图如下图7：图7 N=1024/256周期图谱估计对比从以上结果可以看出N=1024时频谱分辨率明显要高于N=256时的频谱分辨率。（4）、Bartlett平均周期图谱估计当数据长度N=1024及256时，分段数L分别为2和8时信号的功率谱为便于对比，将结果表示如下图8一幅图中：图8 Bartlett平均周期图谱估计结果分析：1、首先数据的长度对分辨率有影响，数据长度N=1024时的频谱分辨率比N=256时的频谱分辨率高；2、分段数L对频谱的分辨率和平滑性（方差）也有很大影响。当数据数目N一定时，L加大，每一段的数据量M就会减小，因此估计方差减小，曲线越平滑，但此时偏移加大，分辨率降低

7、，即估计量的方差和分辨率是一对矛盾，二者的效果可以通过合理选取L互换。3、收获体会1、通过实验，对经典谱估计法有了更深刻的理解，根据课堂中经典谱估计的理论，掌握了周期图法，分段周期图法的具体实现，更加深刻地体会到了理论的原理以及这些估计方法的优缺点，对这些估计方法获得了真正的理解；2、很小的细节也要注意，也需要认真思考，如信号x(n)的产生，通过两种方法的对比，对卷积有了更深入的认识，通过不同方法得到相同结果，学会了从不同的角度分析问题；3、对matlab掌握更好，如自己利用子函数调用的方式实现了对不同分段平均周期图法的实现。函数传递参数分段数L即数据x(n)，子函数代码如下：function

8、 Px_L=Px_mean(L,x)N=length(x); %计算数据长度M=N/L; %每段有M个数据tmp=zeros(1,M); %初始化中间变量，用于保存每段观测数据的功率 for i=1:L %观每段有M=N/L个数据 for j=1:M y(i,j)=x(j+i*M-M); %取出L段数据 end xk=fft(y(i,:),M); %fft变换 Ix=(abs(xk).2)./M; %计算每段的功率 Px=tmp+Ix; %累加每段的功率 end Px_L=Px./L; %平均功率谱二、实验题目二1、问题分析（1）、信号产生首先根据题目给定的信噪比计算正弦信号的幅度，产生所需的

9、正弦信号后与噪声信号叠加形成观测信号。若正弦信号幅度为A，则其功率为：；信噪比：，其中P为噪声功率；由以上公式便可计算出两个正弦信号的幅度。（2）利用MATLAB功率谱分析自相关法首先由观测数据估计出其自相关函数，再解该方程得到模型参数，进而求得信号的功率谱，一般利用Levenson-Durbin递推法求解。matlab中的自相关法谱估计如下： pyulear(xn,p,nfft,fs); 其中， xn是输入观测信号，p是阶数，nfft是观测数据两倍长度的频率采样点，fs是采样频率。Burg法Burg法直接利用观测数据计算AR模型参数，适用于观测数据长度较短的时候。matlab中Burg法谱估

10、计如下： pburg(xn,p,nfft,fs);其中，px1是返回的功率谱，f1是返回的频率分布；xn是输入观测信号，p是阶数，nfft是观测数据两倍长度的频率采样点，fs是采样频率。协方差法协方差法利用到的数据完全一致，相对自相关法没有数据为0的假设。matlab中协方差法谱估计如下： pcov(xn,p,nfft,fs);其中，px1是返回的功率谱，f1是返回的频率分布；xn是输入观测信号，p是阶数，nfft是观测数据两倍长度的频率采样点，fs是采样频率。修正协方差法修正协方差法使用到了前向和后向预测误差对模型参数进行估计。matlab中的修正协方差法谱估计如下： pmcov(xn,p,

11、nfft,fs);其中，px1是返回的功率谱，f1是返回的频率分布；xn是输入观测信号，p是阶数，nfft是观测数据两倍长度的频率采样点，fs是采样频率。2、实验结果与分析（1）数据长度的影响：将四种方法产生的功率谱显示在了一幅图中。L=1024时，如下图9所示。L=256时，结果如下图10所示。 图9 L=1024图10 L=256分析：对比两幅图可以发现：随着数据长度的减小功率谱的分辨率在逐渐降低。即数据长度越大，功率谱分辨率越高，数据长度越小，功率谱分辨率越低。（2）模型阶次的影响：当阶次发生变化时，如p=100时结果如图11所示，p=30时，结果如图12所示。图11 阶次p=100图1

12、2 阶次p=30分析：对比阶次p=100与阶次p=30时的结果可以看出，当阶次太小时，分辨率比较低，无法估计出小信噪比的信号，即小信噪比信号无法在图中观察出来。（3）信噪比的影响：将110Hz正弦信号信噪比降低到10，结果如下图13所示：图13 SNR1=10 SNR2=10对比发现：变化最明显的是自相关法和Burg法，当两个正弦信号信噪比相差不多时，自相关法的估计效果还是较好的；而Burg法效果很差，每次实验的结果不同，经常分辨不出一个正弦频率，出现了谱峰偏移和谱线分裂的现象。而对于自相关法，并不适合于分辨信噪比相差较多的两个正弦信号混合的情况，结果会谱峰偏移，出现观测不到一个正弦信号的频率

13、。如修改信噪比SNR1=30 SNR2=10的情况，结果如下图14所示，仍然分辨不出信噪比较低的信号，即110Hz ，SNR2=10的正弦信号。图13 SNR1=30 SNR2=103、收获体会1、通过实验，加深了对AR模型谱估计的流程，掌握了结合matlab利用自相关法、Burg法、协方差法、修正协方差法估计功率谱，体会了各种方法的性能，优缺点等；2、能更加熟练地通过matlab帮助实现对函数的查找，使用；3、由于考试最近考试较多，时间有限，没有自己编程实现各种现代功率谱估计的方法，在时间空闲时，也会结合课堂讲授的流程图对各种方法尝试自行编程实现。附 matlab源代码实验题目一clc; %

14、清除命令窗口clear; %清空变量N=1024; %数据长度df=1; %设定频率分辨率fs=(N-1)*df; %计算采样率t=0:1/fs:(N-1)/fs; %截取信号的时间段f=0:df:fs; %功率谱估计的频率分辨率和范围w=sqrt(1)*randn(1,N);% 平稳随机信号的获取% 准确获得观测信号是很关键的% 处理方法一：迭代的方式x=w(1) zeros(1,N-1);for i=2:N x(i)=0.8*x(i-1)+w(i);end% 处理方法二：卷积的方式% n=1:N;% h(n)=(0.8).(n-1); %要注意n-1不是n，因为冲击响应是从0开始的% y=

15、conv(w(n),h(n) %共2*N-1项，取前N项即可% 令N=5时，实验结果，取卷积的前N项，可以看出两种方式结果是相同的% x =% 0.6232 1.2976 1.9790 0.5911 0.6849% y =% 0.6232 1.2976 1.9790 0.5911 0.6849 0.3437 0.0131% -0.2978 0.0868% 当数据长度为256时，取前256个数据为x2N2=256;x2=x(1:256);% 理想功率谱计算n=1:N;h(n)=(0.8).(n-1); %要注意n-1不是n，因为冲击响应是从0开始的H=fft(h,N); %对冲击响应做傅里叶变换

16、Px_ideal=(abs(H).2;Px_ideal_db=10*log10(Px_ideal); %以分贝为单位表征% 周期图法估计功率谱% 数据长度为1024时xk=fft(x,N);%傅里叶变换Px_zqt=(abs(xk).2)./N);Px_zqt_db=10*log10(Px_zqt);% 数据长度为256时xk2=fft(x2,N);%傅里叶变换Px_zqt2=(abs(xk2).2)./N2);Px_zqt_db2=10*log10(Px_zqt2);% Px_mean(L,x):平均周期图法,分段数目L，数据为xPx_1024_2=Px_mean(2,x); %参数L=2

17、x(1024个数据)时Px_1024_8=Px_mean(8,x); %参数L=8 x(1024个数据)时Px_256_2=Px_mean(2,x2); %参数L=2 x(256个数据)时Px_256_8=Px_mean(8,x2); %参数L=8 x(256个数据)时Px_1024_2_db=10*log10(Px_mean(2,x); %参数L=2 x(1024个数据)时Px_1024_8_db=10*log10(Px_mean(8,x); %参数L=8 x(1024个数据)时Px_256_2_db=10*log10(Px_mean(2,x2); %参数L=2 x(256个数据)时Px_2

18、56_8_db=10*log10(Px_mean(8,x2); %参数L=8 x(256个数据)时% 显示噪声w(n)与平稳信号x(n)figure;plot(n,w(n); xlabel('n'); ylabel('w(n)');legend('随机噪声w(n)')figure;plot(n,x(n); xlabel('n'); ylabel('x(n)');legend('平稳随机信号x(n)')% 显示理想功率谱figure;plot(f,Px_ideal);xlabel('频率f（H

19、z）');ylabel('功率谱Px_ideal');legend('理想功率谱Px_ideal');figure;plot(f,Px_ideal_db);xlabel('频率f（Hz）');ylabel('功率谱Px_ideal_db（dB）');legend('理想功率谱Px_ideal_db');% 周期图法估计功率谱结果% 数据长度为1024时figure;subplot(2,2,1);plot(f,Px_zqt);xlabel('频率f（Hz）');ylabel('功率谱P

20、x');legend('1024');subplot(2,2,2);plot(f,Px_zqt_db);xlabel('频率f（Hz）');ylabel('功率谱Px（dB）');legend('1024');% 数据长度为256时subplot(2,2,3);plot(f,Px_zqt2);xlabel('频率f（Hz）');ylabel('功率谱Px');legend('256');subplot(2,2,4);plot(f,Px_zqt_db2);xlabel('

21、;频率f（Hz）');ylabel('功率谱Px（dB）');legend('256');% 平均周期图法结果figure;subplot(2,2,1);plot(1:length(Px_1024_2),Px_1024_2);xlabel('频率f（Hz）');ylabel('功率谱Px');legend('L=2 N=1024');subplot(2,2,2);plot(1:length(Px_1024_8),Px_1024_8);xlabel('频率f（Hz）');ylabel('

22、;功率谱Px');legend('L=8 N=1024');subplot(2,2,3);plot(1:length(Px_256_2),Px_256_2);xlabel('频率f（Hz）');ylabel('功率谱Px');legend('L=2 N=256');subplot(2,2,4);plot(1:length(Px_256_8),Px_256_8);xlabel('频率f（Hz）');ylabel('功率谱Px');legend('L=8 N=256');子程序P

23、x_mean.mfunction Px_L=Px_mean(L,x)N=length(x); %计算数据长度M=N/L; %每段有M个数据tmp=zeros(1,M); %初始化中间变量，用于保存每段观测数据的功率 for i=1:L %观每段有M=N/L个数据 for j=1:M y(i,j)=x(j+i*M-M); %取出L段数据 end xk=fft(y(i,:),M); %fft变换 Ix=(abs(xk).2)./M; %计算每段的功率 Px=tmp+Ix; %累加每段的功率 end Px_L=Px./L; %平均功率谱实验题目二clc;clear; fs=1000; % 采样频率 n=0:1/fs:2; %取样数据点 snr1=30; %正弦波d的信噪比snr2=10; %正弦波d2的信噪比dB A1=sqrt(2*10(snr1/10); %正弦波x1的幅度 A2=sqrt(2*10(snr2/10); %正弦波x2的幅度 An=sqrt(1); %噪声w(n)的幅度 x1=A1*sin(2*pi*100*n); %正弦信号x1 x2=A2*sin(2*pi*110*n); %正弦信号x2 w=An*randn(size(n); %噪声信号 w(n) xn=x1+x2+w;