p_classical.m

经典的功率谱估计

clear;close all;clc;
L=1024;%数据总长度为1024
wi=0:2/L:2*(L-1)/L;
randn('state',0);
w=randn(1,L);%产生噪声序列
x=zeros(1,L);%初始化观测序列
x(1)=w(1);
for n=2:L   %产生观测数据序列
    x(n)=0.8*x(n-1)+w(n);
end

a=[1,-0.8];b=1;% 设置原始数据产生的滤波器参数
[h,w1]=freqz(b,a,'whole');% 求出滤波器的频率响应特性
H=(abs(h)).^2;% 取频率响应模的平方与噪声方差（这里为1）相乘即为信号的功率谱密度
figure(1);subplot(2,2,1);
plot(w1/pi,H);grid on;xlim([-0.1,2.1]);
title('理想功率谱密度');
xlabel('角频率ω/π');
ylabel('功率谱密度');
p1=period(x,L,1,L);%调用平均周期图法的函数，分段数为1,相当于没有进行分段平均。因为周期图法是分段数为1的平均周期图法。
subplot(2,2,2);
%p1=p1/min(p1);%找出谱密度中的最小值，作归一化处理
%p1=10*log10(p1);%产生谱密度的dB数据
plot(wi,p1);grid on;xlim([-0.1,2.1]);
title('周期图法（即为平均周期图法分1段的情况）');
xlabel('角频率ω/π');
ylabel('功率谱密度');
p2=period(x,L,2,L);%调用平均周期图法的函数，分段数为2
%p2=p2/min(p2);%找出谱密度中的最小值，作归一化处理
%p2=10*log10(p2);%产生谱密度的dB数据
subplot(2,2,3);
plot(wi,p2);grid on;xlim([-0.1,2.1]);
title('平均周期图法分2段');
xlabel('角频率ω/π');
ylabel('功率谱密度');
p3=period(x,L,8,L);%调用平均周期图法的函数，分段数为4
subplot(2,2,4);
plot(wi,p3);grid on;xlim([-0.1,2.1]);
title('平均周期图法分8段');
xlabel('角频率ω/π');
ylabel('功率谱密度');