功率谱密度估计方法的MATLAB实现

四、 Levinson递推法谱估计程序

1、 源程序

fs=1; %设采样频率为1

N=1000; %数据长度 改变数据长度会导致分辨率的变化； f1=0.2*fs; %第一个sin信号的频率,f1/fs=0.2

f2=0.3*fs; %第二个sin信号的频率，f1/fs=0.2或者0.3

M=16; %滤波器阶数的最大取值,超过则认为代价太大而放弃 L=2*N; %有限长序列进行离散傅里叶变换前，序列补零的长度 n=1:N;

s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s为原始信号

x=awgn(s,10);%x为观测信号，即对原始信号加入白噪声，信噪比10dB figure(1); %画出原始信号和观测信号 subplot(2,1,1);

plot(s,'b'),axis([0 100 -3 3]),xlabel('时间'),ylabel('幅度'),title('原始信号s'); grid; subplot(2,1,2);

plot(x,'r'),axis([0 100 -3 3]),xlabel('时间'),ylabel('幅度'),title('观测信号x'); grid;

%计算自相关函数

rxx = xcorr(x,x,M,'biased');%计算有偏估计自相关函数，长度为-M到M， %共2M+1

r0 = rxx(M+1); %r0为零点上的自相关函数，相对于-M，第M+1个点为零点 R = rxx(M+2:2*M+1);% R为从1到第M个点的自相关函数矩阵 %确定矩阵大小 a = zeros(M,M);

FPE = zeros(1,M);%FPE：最终预测误差，用来估计模型的阶次 var = zeros(1,M); %求初值

a(1,1) = -R(1)/r0;%一阶模型参数 var(1) = (1-(abs(a(1,1)))^2)*r0;%一阶方差

6 / 12

FPE(1) = var(1)*(M+2)/(M); %递推

for p=2:M sum=0;

for k=1:p-1%求a(p,p)

sum=sum+a(p-1,k)*R(p-k); end

a(p,p)=-(R(p)+sum)/var(p-1); for k=1:p-1 %求a(p,k)

a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k); end

var(p)=(1-a(p,p)^2)*var(p-1); %求方差

FPE(p)=var(p)*(M+1+p)/(M+1-p);%求最终预测误差 end

%确定AR模型的最佳阶数

min=FPE(1); %求出FPE最小时对应的阶数 p = 1; for k=2:M if FPE(k)

%功率谱估计

W=0.01:0.01:pi; %功率谱以2*pi为周期，又信号为实信号，只需输出0到PI即可； he=ones(1,length(W)); %length()求向量的长度 for k=1:p

7 / 12

he=he+(a(p,k).*exp(-j*k*W)); end

Pxx=var(p)./((abs(he)).^2); %功率谱函数；

F=W*fs/(pi*2); %将角频率坐标换算成HZ坐标，便于观察；重要！ figure;

plot(F,abs(Pxx))

xlabel('频率/Hz'),ylabel('功率谱P'),title([' AR模型的最佳阶数p=' int2str(p)] ); grid;

2、 仿真结果

8 / 12

五、 Burg法谱估计程序

1、 源程序 fs=1;%设采样频率为1

N=900;%数据长度 改变数据长度会导致分辨率的变化； f1=0.2*fs;%第一个sin信号的频率,f1/fs=0.2

f2=0.3*fs;%第二个sin信号的频率，f1/fs=0.2或者0.3 M=512;%滤波器阶数的最大取值,超过则认为代价太大而放弃 n=1:N;

s = sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s为原始信号

x = awgn(s,10);%x为观测信号，即对原始信号加入白噪声，信噪比10dB for i=1:N ef(1,i)=x(i); eb(1,i)=x(i); end sum=0; for i=1:N

sum=sum+x(i)*x(i);

9 / 12

end

r(1)=sum/N; % Burg递推 for p=2:M

% 求解第p个反射系数 sum1=0; for n=p:N

sum1=sum1+ef(p-1,n)*eb(p-1,n-1); end

sum1=-2*sum1; sum2=0; for n=p:N

sum2=sum2+ef(p-1,n)*ef(p-1,n)+eb(p-1,n-1)*eb(p-1,n-1); end

k(p-1)=sum1/sum2; % 求解预测误差平均功率

r(p)=(1-k(p-1)*k(p-1))*r(p-1); % 求解p阶白噪声方差 q(p)=r(p); % 系数a if p>2

for i=1:p-2

a(p-1,i)=a(p-2,i)+k(p-1)*a(p-2,p-1-i); end end

a(p-1,p-1)=k(p-1); % 求解前向预测误差 for n=p+1:N

ef(p,n)=ef(p-1,n)+k(p-1)*eb(p-1,n-1);

10 / 12