wiener.m

wiener滤波器的matlab程序

clc;
clear all;
a1=-1.352;
a2=1.338;
a3=-0.662;
a4=0.240;
M=200;                  %信号的长度
N=100;                  %滤波器的冲激响应系数的个数
w=randn(1,M);
x=zeros(1,M);           %期望响应
xx=zeros(1,M+4);
xx(1,1)=0;xx(1,2)=0;xx(1,3)=0;xx(1,4)=0;
for i=5:M+4
    xx(1,i)=w(1,i-4)-a1*xx(1,i-1)-a2*xx(1,i-2)-a3*xx(1,i-3)-a4*xx(1,i-4);
end
x(1,1:M)=xx(1,5:M+4);
v=randn(1,M);            %白噪声 
y=x+v;                   %滤波器输入
Re=zeros(1,M);
for i=1:N
    for j=1:M
        if j+i-1<M
            Re(1,i)=y(1,j)*y(1,j+i-1)'/(M-i+1)+Re(1,i);
        end
    end
end
R=zeros(N,N);            %自相关矩阵
for i=1:N
    for j=(i+1):N
        R(i,i)=Re(1,1);
        R(i,j)=Re(1,j-i+1);
        R(j,i)=conj(R(i,j));
    end
end
r=zeros(1,N);             %互相关矩阵
for i=1:N
    for j=1:M
        if j-i+1>0
            r(1,i)=y(1,j)*x(1,j-i+1)'/(M-i+1)+r(1,i);
        end
    end
end
Wopt=inv(R)*r';           %最优抽头权向量
X=zeros(1,M);             %滤波器输出
for n1=1:M
    if n1<N
            for i=0:n1-1
                X(1,n1)=Wopt(i+1,1)*y(1,n1-i)+X(1,n1);
            end
    else
            for i=0:N-1
                X(1,n1)=Wopt(i+1,1)*y(n1-i)+X(1,n1);
            end
    end
end
plot(x);
hold on
plot(X,'r');
hold off