matlab经典、现代功率谱估计

上机作业：1、假设一平稳随机信号为()()()0.81x n x n w n =−+，其中 是均值为0，方差为1的白噪声，数据长度为1024。

（1）、产生符合要求的)(n w 和)(n x ；（2）、给出信号)(n x 的理想功率谱；（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(e jw )=σw 2|H(e jw )|2。

周期图方法：直接对观测数据做FFT 变换，变换的结果取模的平方再除以数据长度，作为估计的功率谱。

256个观测点时可以对原观测数据以4为间隔提取得到。

Bartlett 法：将L 组独立的观测数据分别求周期图，再将L 个周期图求平均作为信号的功率谱估计。

L 组数据可以通过对原观测数据以L 为间隔提取得到。

二、MATLAB 实现程序及注解 clc;clear;close all;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)for i=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));grid on;title('理想功率谱');xlabel('频率'); ylabel('功率db');%% 周期图法%1024个观测点Pxx=abs(fft(x)).^2/N; %周期图公式Pxx=10*log10(Pxx(index+1)); %化为dbfigure;plot(k,Pxx);grid on;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);grid on;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)));grid on;title('N=1024,L=2');xlabel('频率'); ylabel('功率db');%分段L=8L1=8;x3=zeros(L1,N/L1); %产生L1行,N/L1列的矩阵用以存储分组的数据for i=1:L1x3(i,:)=x(i:L1:N); %将原始数据分为8组endPxx3=zeros(L1,1024); %产生L1行,1024列矩阵用以存储分组的周期图for i=1:L1Pxx3(i,:)=abs(fft(x3(i,:),1024)).^2/length(x3(i,:)); %分别求周期图,结果保存在Pxx3中,FFT长度为1024endfor i=1:1024Pxx3_m(i)=sum(Pxx3(:,i))/L1; %求平均endsubplot(2,2,2),plot(k,10*log10(Pxx3_m(index+1)));grid on;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)));grid on;title('N=256,L=2');xlabel('频率'); ylabel('功率db');%分段L=8L3=8;x4=zeros(L3,length(x)/L3);for i=1:L3x4(i,:)=x(i:L3:length(x)); %将原始数据分为8组endPxx4=zeros(L3,1024);for i=1:L3Pxx4(i,:)=abs(fft(x4(i,:),1024)).^2/length(x4(i,:)); %分别求周期图,FFT长度为1024endfor i=1:1024Pxx4_m(i)=sum(Pxx4(:,i))/L3; %求平均endsubplot(2,2,4),plot(k,10*log10(Pxx4_m(index+1)));grid on;title('N=256,L=8');xlabel('频率'); ylabel('功率db');三、结果及分析图1 理想功率谱图2 周期图1024点及256点从上图可以看出，周期图法得到的功率谱估计，谱线的起伏较大，即估计所得的均方误差较大。

当N增加时，摆动的频率加快，而摆动的幅度变化不大。

且N=1024时，谱的分辨率较N=256时大。

图3 Bartlett平均周期图法由上图可以看出，采用平均处理后，谱线上下摆动的幅度减小（即均方误差有所降低），曲线的平滑性也较周期图法好。

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，对其在一个周期内的平方进行积分再除以周期长度可得，平均功率P s=A2/2。

白噪声的平均功率为1，所以又SNR的计算式可得题中两个正弦信号的幅值分别为sqrt(20)和sqrt(2000)。

题中，所给观测信号为白噪声混有两个正弦信号，所以功率谱中应有两个明显的尖峰，所以选择AR模型来进行估计。

AR模型的估计方法有自相关法，Burg法，协方差法和修正协方差法等。

本次编程采用matlab自带的功率谱分析函数。

自相关，Burg法，协方差法，修正协方差法分别为pyulear()，pburg()，pcov()，pmcov()。

二、MATlAB实现程序及注解clc;clear;close all;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));grid on;title('自相关法');[Pxx1,f1]=pcov(x,p,nfft,fs);figure,plot(f1,10*log10(Pxx1));grid on;title('协方差法');[Pxx2,f2]=pmcov(x,p,nfft,fs);figure,plot(f2,10*log10(Pxx2));grid on;title('修正协方差法');[Pxx3,f3]=pburg(x,p,nfft,fs);figure,plot(f3,10*log10(Pxx3));grid on;title('Burg法');三、结果及分析图1 50阶AR模型谱估计从上图可以看出，采用参数建模的谱估计方法得到的功率谱曲线平滑（方差小），分辨率高，可以明显地观察到两个谱峰。

图2 20阶AR模型谱估计降低模型阶次后，可以发现，谱的分辨率降低（两个谱峰几乎变成一个谱峰），但是曲线平滑性变好（估计误差降低）。

图3 50阶AR模型（4倍观测数据）谱估计与图1相比，图3几乎没有什么变化，这是因为AR模型隐含着数据的外推。