matlab功率谱估计

功率谱估计及其MATLAB仿真

1经典功率谱估计

经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。

1.1相关函数法(BT法)

该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。

Matlab代码示例1(Btfangfa.M)：

Fs=500;%采样频率

n=0:1/Fs:1;

xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));%产生含有噪声的序列

nfft=512;

cxn=xcorr(xn,'unbiased');%计算序列的自相关函数

CXk=fft(cxn,nfft);

Pxx=abs(CXk);

index=0:round(nfft/2-1); %Round towards nearest integer.

k=index*Fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

figure(1);

plot(k,plot_Pxx);

结果如下：

1.2周期图法(periodogram)

周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换得X(k),然后再取其幅值的平方,并除以N,作为序列x(n)真实功率谱的估计。Matlab代码示例2(PEriod.M)：

Fs=600;

n=0:1/Fs:1;

xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));

window=boxcar(length(xn));%矩形窗

nfft=512;

[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法

figure(1);

plot(f,10*log10(Pxx));

window=boxcar(length(xn));%矩形窗

nfft=1024;

[Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法

figure(2);

plot(f,10*log10(Pxx));

结果如下：

1.3平均周期图法和平滑平均周期图法

对于周期图的功率谱估计,当数据长度N太大时,谱曲线起伏加剧,若N太小,谱的分辨率又不好,因此需要改进。两种改进的估计法是平均周期图法和平滑平均周期图法。Bartlett 法:Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。

Matlab代码示例3(Bartlett.M)：

fs=600;

n=0:1/fs:1;

xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n));

nfft=512;

window=hamming(nfft);%矩形窗

noverlap=0;%数据无重叠

p=0.9;%置信概率

[Pxx,Pxxc]=psd(xn,nfft,fs,window,noverlap,p);

index=0:round(nfft/2-1);

k=index*fs/nfft;

plot_Pxx=10*log10(Pxx(index+1));

plot_Pxxc=10*log10(Pxxc(index+1));

figure(1)

plot(k,plot_Pxx);

figure(2)

plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]);

结果如下：

Welch法:Welch法对Bartlett法进行了两方面的修正,一是选择适当的窗函数w(n),并在周期图计算前直接加进去,加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时,可使各段之间有重叠,这样会使方差减小。

Matlab代码示例4(Welch.M)：

Fs=600;

n=0:1/Fs:1;

xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));

nfft=512;

window=boxcar(100);%矩形窗

window1=hamming(100);%海明窗

window2=blackman(100);%blackman窗

noverlap=20;%数据无重叠

range='half';%频率间隔为[0 Fs/2],计算一半的频率

[Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range);

[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range);

[Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range);

plot_Pxx=10*log10(Pxx);

plot_Pxx1=10*log10(Pxx1);

plot_Pxx2=10*log10(Pxx2);

figure(1)

plot(f,plot_Pxx);

figure(2)

plot(f,plot_Pxx1);

figure(3)

plot(f,plot_Pxx2);

结果如下：