matlab经典、现代功率谱估计

上机作业：1、假设一平稳随机信号为()()()0.81xnxnwn=−+，其中是均值为0，方差为1的白噪声，数据长度为1024。（1）、产生符合要求的)(nw和)(nx；（2）、给出信号)(nx的理想功率谱；（3）、编写周期图谱估计函数，估计数据长度N=1024及256时信号功率谱，分析估计效果。（4）、编写Bartlett平均周期图函数，估计当数据长度N=1024及256时，分段数L分别为2和8时信号的功率谱，分析估计效果。一、解题思路w(n)可以通过随机序列randn(1,N)来产生，x(n)可以通过对w(n)滤波产生（由递推式可得系统的传递函数），也可以直接由递推式迭代产生。由于线性系统的输出功率谱等于输入功率谱乘以传递函数模的平方，X(n)可以看做w(n)通过一线性系统的输出，H(z)=1/(1-0.8z)。所以x(n)的理想功率谱P(ejw)=σw2|H(ejw)|2。周期图方法：直接对观测数据做FFT变换，变换的结果取模的平方再除以数据长度，作为估计的功率谱。256个观测点时可以对原观测数据以4为间隔提取得到。Bartlett法：将L组独立的观测数据分别求周期图，再将L个周期图求平均作为信号的功率谱估计。L组数据可以通过对原观测数据以L为间隔提取得到。二、MATLAB实现程序及注解clc;clear;closeall;Fs=500;%采样率N=1024;%观测数据w=sqrt(1)+randn(1,N);%0均值,方差为1的白噪声,长度1024x=[w(1)zeros(1,N-1)];%初始化x(n),长度1024,x(1)=w(1)fori=2:Nx(i)=0.8*x(i-1)+w(i);%迭代产生观测数据x(n)end%%理想功率谱[h,w1]=freqz(x);figure,plot(w1*500/(2*pi),10*log10(abs(h).^2));gridon;title('理想功率谱');xlabel('频率');ylabel('功率db');%%周期图法%1024个观测点Pxx=abs(fft(x)).^2/N;%周期图公式Pxx=10*log10(Pxx(index+1));%化为dbfigure;plot(k,Pxx);gridon;title('周期图1024点');xlabel('频率');ylabel('功率db');%周期图256个观测点x1=x(1:4:N);Pxx1=abs(fft(x1,1024)).^2/N;Pxx1=10*log10(Pxx1(index+1));%化为dbfigure;plot(k,Pxx1);gridon;title('周期图256点');xlabel('频率');ylabel('功率db');%%Bartlett平均周期图,N=1024%分段L=2L=2;x_21=x(1:L:N);x_22=x(2:L:N);Pxx_21=abs(fft(x_21,1024)).^2/length(x_21);Pxx_22=abs(fft(x_22,1024)).^2/length(x_22);Pxx_2=(Pxx_21+Pxx_22)/L;figure;subplot(2,2,1),plot(k,10*log10(Pxx_2(index+1)));gridon;title('N=1024,L=2');xlabel('频率');ylabel('功率db');%分段L=8L1=8;x3=zeros(L1,N/L1);%产生L1行,N/L1列的矩阵用以存储分组的数据fori=1:L1x3(i,:)=x(i:L1:N);%将原始数据分为8组endPxx3=zeros(L1,1024);%产生L1行,1024列矩阵用以存储分组的周期图fori=1:L1Pxx3(i,:)=abs(fft(x3(i,:),1024)).^2/length(x3(i,:));%分别求周期图,结果保存在Pxx3中,FFT长度为1024endfori=1:1024Pxx3_m(i)=sum(Pxx3(:,i))/L1;%求平均endsubplot(2,2,2),plot(k,10*log10(Pxx3_m(index+1)));gridon;title('N=1024,L=8');xlabel('频率');ylabel('功率db');%%Bartlett平均周期图,N=256,求法同上%分段L=2,分别计算周期图,再取平均x=x(1:4:N);L2=2;x_31=x(1:L2:length(x));x_32=x(2:L2:length(x));Pxx_31=abs(fft(x_31,1024)).^2/length(x_31);Pxx_32=abs(fft(x_32,1024)).^2/length(x_32);Pxx_3=(Pxx_31+Pxx_32)/L2;subplot(2,2,3),plot(k,10*log10(Pxx_3(index+1)));gridon;title('N=256,L=2');xlabel('频率');ylabel('功率db');%分段L=8L3=8;x4=zeros(L3,length(x)/L3);fori=1:L3x4(i,:)=x(i:L3:length(x));%将原始数据分为8组endPxx4=zeros(L3,1024);fori=1:L3Pxx4(i,:)=abs(fft(x4(i,:),1024)).^2/length(x4(i,:));%分别求周期图,FFT长度为1024endfori=1:1024Pxx4_m(i)=sum(Pxx4(:,i))/L3;%求平均endsubplot(2,2,4),plot(k,10*log10(Pxx4_m(index+1)));gridon;title('N=256,L=8');xlabel('频率');ylabel('功率db');三、结果及分析图1理想功率谱图2周期图1024点及256点从上图可以看出，周期图法得到的功率谱估计，谱线的起伏较大，即估计所得的均方误差较大。当N增加时，摆动的频率加快，而摆动的幅度变化不大。且N=1024时，谱的分辨率较N=256时大。图3Bartlett平均周期图法由上图可以看出，采用平均处理后，谱线上下摆动的幅度减小（即均方误差有所降低），曲线的平滑性也较周期图法好。N相同时，分段越多，方差越小，曲线越平滑。这是因为，N一定时，L加大，每一段的数据量就会相应减少，因此估计方差减小，偏移加大，从而分辨率降低。N一定时，L与每段数据量相互矛盾，需择中选取。综上，Bartlett法相对于周期图法来说，较好地减小了估计误差。2、假设均值为0，方差为1的白噪声w(n)中混有两个正弦信号，该正弦信号的频率分别为100Hz和110Hz，信噪比分别为10dB和30dB，初始相位都为0，采样频率为1000Hz。（1）、采用自相关法、Burg法、协方差法、修正协方差法估计功率谱，分析数据长度和模型阶次对估计结果的影响（可采用MATLAB自带的功率谱分析函数）。（2）、调整正弦信号信噪比，分析信噪比的降低对估计效果的影响。一、解题思路信噪比为信号平均功率与噪声的平均功率之比，即SNR=10lgS/N，单位dB。假设正弦信号的幅度为A，对其在一个周期内的平方进行积分再除以周期长度可得，平均功率Ps=A2/2。白噪声的平均功率为1，所以又SNR的计算式可得题中两个正弦信号的幅值分别为sqrt(20)和sqrt(2000)。题中，所给观测信号为白噪声混有两个正弦信号，所以功率谱中应有两个明显的尖峰，所以选择AR模型来进行估计。AR模型的估计方法有自相关法，Burg法，协方差法和修正协方差法等。本次编程采用matlab自带的功率谱分析函数。自相关，Burg法，协方差法，修正协方差法分别为pyulear()，pburg()，pcov()，pmcov()。二、MATlAB实现程序及注解clc;clear;closeall;fs=1000;%采样率1000N=1;%改变数据长度p=50;%AR模型阶数nfft=512;%fft长度t=0:1/fs:N;wn=sqrt(1)+randn(1,N*fs+1);%白噪声,均值0,方差1s1=sqrt(20)*sin(2*pi*100*t);%正弦信号1,信噪比10dbs2=sqrt(2000)*sin(2*pi*110*t);%正弦信号2,信噪比30dbx=s1+s2+wn;%观测数据%figure,plot(t,x);x1=xcorr(x,'biased');[Pxx,f]=pyulear(x1,p,nfft,fs);%Yule-Walker方程figure,plot(f,10*log10(Pxx));gridon;title('自相关法');[Pxx1,f1]=pcov(x,p,nfft,fs);figure,plot(f1,10*log10(Pxx1));gridon;title('协方差法');[Pxx2,f2]=pmcov(x,p,nfft,fs);figure,plot(f2,10*log10(Pxx2));gridon;title('修正协方差法');[Pxx3,f3]=pburg(x,p,nfft,fs);figure,plot(f3,10*log10(Pxx3));gridon;title('Burg法');三、结果及分析图150阶AR模型谱估计从上图可以看出，采用参数建模的谱估计方法得到的功率谱曲线平滑（方差小），分辨率高，可以明显地观察到两个谱峰。图220阶AR模型谱估计降低模型阶次后，可以发现，谱的分辨率降低（两个谱峰几乎变成一个谱峰），但是曲线平滑性变好（估计误差降低）。图350阶AR模型（4倍观测数据）谱估计与图1相比，图3几乎没有什么变化，这是因为AR模型隐含着数据的外推。在经典谱估计中，对于观测数据以外的数据均假设为0（相当于时域加窗），从而造成谱的分辨率降低。AR谱估计中，由于隐含自相关函数的外推，而使分辨率大大提高。从而观测数据的多少对分辨率影响不大，对估计误差会有影响。数据越多，估计误差越小，谱线波动越小，这一点图中可以很明显地看出。图420阶AR模型（增大SNR）谱估计上图是将原信噪比变为50dB和60dB后得到的功率谱估计，与图2相比，图中能够看出两个谱峰，即分辨率较高。所以相同阶次时，信噪比越高，分辨率越高，信噪比降低则会导致分辨率降低。对于已知信噪比的观测信号来说，若信噪比低则选用较高的阶次来得到较好地分辨率；信噪比高，则相应地可以减小阶次，以得到较小的误差。