iir_lms.m

用总体最小二乘算法进行ARMA功率谱估计

clear all;
N=64;
n=1:N;
%%%%%设置观测信号%%%%%%%%
v=randn(1,N);
u(n)=12*sin(2*pi*0.22*n);
x(n)=u(n)+v(n);
%%%%%%%求解陷波器参数%%%%%%%%
s=zeros(3,N);
r=zeros(2,N);
for i=1:N;
    S(1,i)=x(i);
    R(1,i)=x(i);
end
C=zeros(2,N);
D=zeros(2,N);
p=0.9;
K0=zeros(1,N);
%%%%%%LMS算法%%%%%%%%
for i=2:N;
    C(1,i)=p*C(1,i-1)+S(1,i)*R(1,i-1);
    D(1,i)=p*D(1,i-1)+1/2*(S(1,i)*S(1,i)+R(1,i-1)*R(1,i-1));
    K0(i)=-C(1,i)/D(1,i);
end
K1=zeros(1,N);
for i=2:N;
    S(2,i)=S(1,i)+K0(i)*R(1,i-1);
    R(2,i)=K0(i)*S(1,i)+R(1,i-1);
    C(2,i)=p*C(2,i-1)+S(2,i)*R(2,i-1);
    D(2,i)=p*D(2,i-1)+1/2*(S(2,i)*S(2,i)+R(2,i-1)*R(2,i-1));
    K1(i)=-C(2,i)/D(2,i);
    S(3,i)=S(2,i)+K1(i)*R(2,i);
end
q=0.99;
%%%%%%%%画出某时刻陷波器的频率响应%%%%%%%
b=[1,K0(N)*(1+K1(N)),K1(N)];
a=[1,K0(N)*(1+q*K1(N)),q*K1(N)];
freqz(b,a,N);
%%%%%%%用陷波器构成普线增强器%%%%%%%%
v1=S(3,:);
u1=x-v1;
figure(2)
plot(n,u,n,u1,'r');
legend('输入信号u（n)','滤波输出u1(n)');