gonglvpuguji.m

各种功率谱估计算法!包括AR谱估计,BURG算法,YULE-WALK方程

%% file name:数字信号处理2大作业---功率谱估计
% Author: 李军    
% 西安电子科技大学，研2-061班，学号：0610120520，Tel:[phone]

%% 建立二阶AR模型，参数初始化
% x(n)=-a1*x(n-1)-a2*x(n-2)+v(n)
clc
clear all
close all
L=5000;%数据长度
s(1)=2;
s(2)=2;
v(1)=0;
v(2)=0;
gamma=0.35;%噪声方差
v(3:L)=sqrt(gamma)*randn(1,L-2);
for i=3:L
    s(i)=-0.1*s(i-1)+0.8*s(i-2)+v(i);
end
figure,plot(s);
legend('二阶AR模型输入信号');
xlabel('i');
ylabel('信号s(i)');

%% 利用FFT求解功率谱估计
xiangguan=1/L*xcorr(s);%求序列自相关
pw1=fft(xiangguan);
length1=length(pw1);
n1=0:length1-1;
ommiga=2*pi*n1/length1;
figure,plot(ommiga,abs(pw1));
legend('利用FFT相关图法来估计功率谱');
xlabel('角频率\omega');
ylabel('功率谱P_x_x(\omega)');

%% 利用AR模型的Yule-Walker方程求解模型参数{a1,a2}
xiangguan=xcorr(s);
fai=toeplitz(1/L*xiangguan(L:(L+2)));
fai1=fai(2:3,2:3);
fai2=(fai(2:3,1));
a=-inv(fai1)*fai2;%估计模型参数{a1,a2}
disp('利用AR模型的Yule-Walker方程求解模型参数{a1,a2}:');disp(a);
gammaguji=fai(1,1:3)*[1;a];%估计噪声方差
disp('噪声方差：');disp(gammaguji);
x=0:0.1:(2*pi-0.1);
i=1;
Pw=zeros(1,length(x));%计算功率谱
for omiga=0:0.1:(2*pi-0.1)    
    Pw(i)=gammaguji/(abs(1+0.1*exp(-j*omiga*1)+(-0.8)*exp(-j*omiga*2)))^2;
    i=i+1;
end
figure,plot(x,Pw);
title('利用AR模型的Yule-Walker方程的功率谱估计');
legend('输入信号的功率谱估计');
xlabel('角频率\omega');
ylabel('功率谱P_x_x(\omega)');

%% 使用LMS算法来估计线性预测器的模型参数{a1,a2}
% LMS算法：
% W(j+1)=W(j)+2*u*e(j)*X(j)
% y(j)=W(j)'*X(j)
% e(j)=d(j)-y(j)
% 线性预测器: x~(n)=-a1*x(n-1)-a2*x(n-2)
% 线性预测误差：e(n)=x(n)-x~(n)=x(n)-a1*x(n-1)-a2*x(n-2)
% 将x(n)作为输入，e(n)作为输出,则线性预测误差输出应该为白噪声w(n)
w=zeros(L,2);
w(2,1:2)=[1 1];%模型参数赋初值
mu=0.0025;%步长因子
y=zeros(1,L);
e=zeros(1,L);
for j=3:L
    y(j)=s(j-1:-1:j-2)*w(j-1,1:2)';%注意这里的w与模型参数相差‘-’号
    e(j)=s(j)-y(j);
    w(j,1:2)=w(j-1,1:2)+2*mu*e(j)*s(j-1:-1:j-2);
end
w=-w';%refesh
figure,plot(w(1,:));
title('LMS算法模型参数\{a1\}逼近曲线');
xlabel('迭代次数 n');
ylabel('待估计参数 a1');
legend('模型参数\{a1\}迭代轨迹');
figure,plot(w(2,:));
title('LMS算法模型参数\{a2\}逼近曲线');
xlabel('迭代次数 n');
ylabel('待估计参数 a2');
legend('模型参数\{a2\}迭代轨迹');
figure,plot(w(1,:),w(2,:));
title('LMS算法模型参数\{a1,a2\}逼近曲线');
xlabel('待估计参数 a1');
ylabel('待估计参数 a2');
legend('模型参数\{a1,a2\}迭代轨迹');
figure,plot(v,'r');
hold on
plot(e,'g--');
legend('原始白噪声','估计白噪声');
x=0:0.1:(2*pi-0.1);
i=1;
jjj=sqrt(-1);
Pw0=zeros(1,length(x));%计算功率谱
for omiga=0:0.1:(2*pi-0.1)    
    Pw0(i)=gamma/(abs(1+w(1,L)*exp(-jjj*omiga*1)+w(2,L)*exp(-jjj*omiga*2)))^2;
    i=i+1;
end
figure,plot(x,Pw0);
title('LMS算法估计功率谱');
legend('使用LMS算法来估计功率谱');
xlabel('角频率\omega');
ylabel('功率谱P_x_x(\omega)');
disp('使用LMS算法来估计线性预测器的模型参数{a1,a2}：');disp(w(:,L));

%% 使用RLS算法来估计线性预测器的模型参数{a1,a2}
% RLS算法简介：
% y(k)=X(k)*W';
% e(k)=d(k)-y(k)
% inv(R(k))=inv(R(k-1))-inv(R(k-1))*X(k)'*X(k)*inv(R(k-1))/(1+X(k)*inv(R(k-1))*X(k)')
% W(k)'=W(k-1)'+inv(R(k-1))*X(k)'*e(k)
w2=zeros(2,L);
w2(1:2,2)=[1;1];
R=eye(2,2);
y2=zeros(1,L);
e2=zeros(1,L);
for j=3:L
    y2(j)=s(j-1:-1:j-2)*w2(:,j-1);%注意这里的w2与模型参数相差‘-’号
    R=R-R*(s(j-1:-1:j-2))'*s(j-1:-1:j-2)*R/(1+s(j-1:-1:j-2)*R*(s(j-1:-1:j-2))');
    e2(j)=s(j)-y2(j);
    w2(:,j)=w2(:,j-1)+R*(s(j-1:-1:j-2))'*e2(j);
end
w2=-w2;
figure,plot(v,'g--');
hold on
plot(e2,'b-.')
legend('原始白噪声','估计白噪声');
x=0:0.1:(2*pi-0.1);
i=1;
Pw1=zeros(1,length(x));%计算功率谱
for omiga=0:0.1:(2*pi-0.1)    
    Pw1(i)=gamma/(abs(1+w2(1,L)*exp(-jjj*omiga*1)+w2(2,L)*exp(-jjj*omiga*2)))^2;
    i=i+1;
end
figure,plot(x,Pw1);
title('RLS算法估计功率谱');
legend('使用RLS算法来估计功率谱');
xlabel('角频率\omega');
ylabel('功率谱P_x_x(\omega)');
disp('使用RLS算法来估计线性预测器的模型参数{a1,a2}：');disp(w2(:,L));
figure,plot(w2(1,:));
title('RLS算法模型参数\{a1\}逼近曲线');
xlabel('迭代次数 n');
ylabel('待估计参数 a1');
legend('模型参数\{a1\}迭代轨迹');
figure,plot(w2(2,:));
title('RLS算法模型参数\{a2\}逼近曲线');
xlabel('迭代次数 n');
ylabel('待估计参数 a2');
legend('模型参数\{a2\}迭代轨迹');
figure,plot(w2(1,:),w2(2,:));
title('RLS算法模型参数\{a1,a2\}逼近曲线');
xlabel('待估计参数 a1');
ylabel('待估计参数 a2');
legend('模型参数\{a1,a2\}迭代轨迹');

%% 使用Levinson-Durbin算法来估计线性预测器的模型参数{a1,a2}
xiangguan=1/L*xcorr(s);
a11=-xiangguan(L+1)/xiangguan(L);
sigma1=(1-(abs(a11)).^2)*xiangguan(L);
a22=-(xiangguan(L+2)+a11*xiangguan(L+1))/sigma1;
a21=a11+a22*a11;
sigma2=(1-(abs(a22)).^2)*sigma1;
w3=[a21;a22];
x=0:0.1:(2*pi-0.1);
i=1;
Pw2=zeros(1,length(x));%计算功率谱
for omiga=0:0.1:(2*pi-0.1)    
    Pw2(i)=sigma2/(abs(1+a21*exp(-jjj*omiga*1)+a22*exp(-jjj*omiga*2)))^2;
    i=i+1;
end
figure,plot(x,Pw2);
title('Levinson-Durbin算法估计功率谱');
legend('使用Levinson-Durbin算法来估计功率谱');
xlabel('角频率\omega');
ylabel('功率谱P_x_x(\omega)');
disp('使用Levinson-Durbin算法来估计线性预测器的模型参数{a1,a2}：');disp(w3);

a=-xiangguan(L+1)/xiangguan(L);
deta=(1-abs(a).^2)*xiangguan(L);
N=4;%阶数长度
figure;
for k=2:N
b=a;
a=ones(k,1);
a(k,1)=-(xiangguan(L+k)+xiangguan(L+k-1:-1:L+1)*b)/deta;
for j=1:k-1
a(j,1)=b(j,1)+a(k,1)*b(k-j,1);
end
deta=(1-abs(a(k,1)).^2)/deta;
x=0:0.1:(2*pi-0.1);
i=1;
Pw2=zeros(1,length(x));%计算功率谱
for omiga=0:0.1:(2*pi-0.1)    
    Pw22(i)=deta/(abs(1+sum(exp(-jjj*omiga*(1:length(a)))*a)))^2;
    i=i+1;
end
plot(x,Pw22);
hold on
end
title('阶数长度变化时对功率谱的影响');
xlabel('角频率\omega');
ylabel('功率谱P_x_x(\omega)');

%% 使用Burg算法来估计线性预测器的模型参数{a1,a2}
e0=s;
b0=s;
sigma02=1/L*sum(s.^2);
k1=-(2*sum(s(2:L)+s(1:L-1)))/sum(s(2:L).^2+s(1:L-1).^2);
e1(1)=x(1);
e1(2:L)=s(2:L)+k1.*s(1:L-1);
b1(1)=k1.*x(1);
b1(2:L)=s(1:L-1)+k1.*s(2:L);
sigma12=(1-k1^2)*sigma02;
k2=-2*(sum(e1(2:L).*b1(1:L-1)))/sum(e1(2:L).^2+b1(1:L-1).^2);
e2(1)=x(1);
e2(2:L)=e1(2:L)+k2.*b1(1:L-1);
b2(1)=k1.*x(1);
b2(2:L)=b1(1:L-1)+k2.*e1(2:L);
sigma22=(1-k2.^2)*sigma12;
x=0:0.1:(2*pi-0.1);
i=1;
Pw3=zeros(1,length(x));%计算功率谱
for omiga=0:0.1:(2*pi-0.1)    
    Pw3(i)=sigma22/(abs(1+a21*exp(-jjj*omiga*1)+a22*exp(-jjj*omiga*2)))^2;
    i=i+1;
end
figure,plot(x,Pw3);
title('Burg算法估计功率谱');
legend('使用Burg算法来估计功率谱');
xlabel('角频率\omega');
ylabel('功率谱P_x_x(\omega)');
disp('使用Burg算法来估计线性预测器的模型参数{a1,a2}：');disp([a21;a22]);
disp('使用Burg算法来估计线性预测器的白噪声功率谱密度估计：');disp(sigma22);
psd