p_per.m

谱估计里的周期图仿真程序

%%%
clear all;close all;clc;
N=32;
M=10;
a1=-0.850848;
sigma2=0.101043;
f1=0.05;
f2=0.40;
f3=0.42;
f=-0.50:0.0001:0.50;
%caculate x_n
u_R= normrnd(0,sigma2/2,1,N);
u_I= normrnd(0,sigma2/2,1,N);
u_n=u_R+j*u_I;
for n=1:N-1
    z(1)=u_n(1);
    z(n+1)=-a1*z(n)+u_n(n+1);
end
for n=1:N-1
    x(1)=6+z(1);
    x(n+1)=2*cos(2*pi*f1*n)+2*cos(2*pi*f2*n)+cos(2*pi*f3*n)+z(n+1);
end
%周期图谱估计法
s=0;
for n=1:N
    s=s+x(n).*exp(-j.*2.*pi.*f.*n);
end
P_per=abs(s).^2./N;
P_per_dB=10.*log10(P_per);
%caculate w_b(k)
k=1:N;
w_B=1-abs(k)./N;
%caculate r_xx 有偏自相关
for k=1:N
    s=0;
    for n=1:N-k
        s=s+(x(n+k)*conj(x(n)));
    end
    r_xx(k)=1/N*s;
end
%布莱克曼-杜基法
P1=0;
P2=0;
for k=1:M+1
    P1=P1+w_B(k).*r_xx(k).*exp(-j.*2.*pi.*f.*k);
end
for k=-M-1:-1
    P2=P2+w_B(-k).*conj(r_xx(-k)).*exp(-j.*2.*pi.*f.*k);
end
P_bt=P1+P2;
P_bt_dB=10.*log10(P_bt);
%画图
figure
plot(f,P_per_dB);
axis([-0.50,0.50,-30.00,30.00]);
xlabel('频率');
ylabel('PSD(dB)');
title('周期图');
figure
plot(f,P_bt_dB);
axis([-0.50,0.50,-30.00,30.00]);
xlabel('频率');
ylabel('PSD(dB)');
title('布莱克曼-杜基');