power_est.m

在Matlab环境下实现的功率谱估计程序

clear;
%调用数据文件test.mat
load test x;
N=4096;
Fn=-0.5:1/N:0.5-1/N;

%直接法
n=length(x);
X=fft(x,N);
P_per=(abs(X)).^2/n;
%归一化频谱
P_per=abs(P_per+eps)/abs(max(P_per));
%绘出估计的功率谱
subplot(2,2,1);
f1=fftshift(10*log10(P_per));
plot(Fn,f1);
axis([-0.5 0.5 -50 0]);grid on;
title('直接法功率谱');

% 间接法
%先求得自相关函数
%然后由自相关函数的傅里叶变换得到功率谱
M=32;
%求自相关函数
rx=xcorr(x,M,'biased');
P_bt=fft(rx,N);
P_bt=abs(P_bt+eps)/abs(max(P_bt));
%绘出估计的功率谱
subplot(2,2,2);
f2=fftshift(10*log10(P_bt));
plot(Fn,f2);
axis([-0.5 0.5 -50 0]);grid on;
title('间接法功率谱');


% Bartlett法
Mb=32;
S=size(x);
L=floor(S(2)/Mb);
%构造矩形窗
w=ones(1,Mb);
PB_per0=zeros(1,N);
for i=1:L
    xnew=x((i-1)*Mb+1:i*Mb);  %分成L段，每段Mb个
    xnew=xnew.*w;
    X=fft(xnew,N);
    PB_per=(abs(X)).^2/Mb;
    PB_per0=(PB_per0+PB_per)/L;
end
PB_per0=abs(PB_per0+eps)/abs(max(PB_per0));
f3=fftshift(10*log10(PB_per0));
subplot(2,2,3);
plot(Fn,f3);
axis([-0.5 0.5 -50 0]);grid on;
title('Bartlett法功率谱');


% Welch法
Mw=32;
S=size(x);
L=floor((S(2)-Mw/2)/Mw*2);
%构造汉明窗
n=1:Mw;
w_hamming=0.54-0.46*cos(2*pi*(n-1)/N);
PB_per1=zeros(1,N);
for i=1:L
    xnew=x((i-1)*Mw/2+1:(i+1)*Mw/2);  %分成L段，每段Mw个，重叠Mw/2个
    xnew=xnew.*w_hamming;
    X=fft(xnew,N);
    PB_per=(abs(X)).^2/sum(w_hamming.*w_hamming);
    PB_per1=(PB_per1+PB_per)/L;
end
PB_per1=abs(PB_per1+eps)/abs(max(PB_per1));
f4=fftshift(10*log10(PB_per1));
subplot(2,2,4);
plot(Fn,f4);
axis([-0.5 0.5 -50 0]);grid on;
title('Welch法功率谱');