wienerfilter2.m

基于维纳滤波器算法,对加噪信号进行滤波.

%维纳滤波程序^_^  李军 2006.10.23
clear all
close all
clc
%选择输入点个数和滤波器阶数----------------------------------------------
L=200;
N=20;
% L=input('输入采样点个数 L= ');
% N=input('输入FIR滤波器阶数 N= ');

%原始信号s(n)，输入信号x(n)=s(n)+v(n)，噪声v(n),信号估计sguji(n)-------
w=1*randn(1,L);
s(1)=w(1);
for n=2:L
    s(n)=s(n-1)+w(n);
end

v=1*randn(1,L);
x=s+v;
figure,plot(x,'r-');
hold on
plot(s,'b');
hold on

%求x(n)的自相关序列和x(n)与s(n)的互相关序列------------------------------
faixx=zeros(1,L);
faixs=zeros(1,L);
for m=1:N
    for n=1:L-m+1
        faixx(1,m)=faixx(1,m)+x(n)*x(n+m-1)/(L-m+1);
        faixs(1,m)=faixs(1,m)+x(n)*s(n+m-1)/(L-m+1);
    end
end

% b=faixx;
% a=faixx+faivv;
% [r,p,c]=residuez(b,a);

%因为是平稳随机序列，则可以将自相关序列扩展为自相关矩阵------------------
%求出h(n)
fai_xx=toeplitz(faixx);
hn=inv(fai_xx)*faixs';

%由h(n)求出对原始信号s(n) 的估计 sguji(n)
sguji=zeros(1,L);
for i=1:L
    if (i<N|i==N)
        for j=1:i
        sguji(1,i)=sguji(1,i)+hn(j,1)*x(1,i-j+1);
        end
    else
        for j=1:N
        sguji(1,i)=sguji(1,i)+hn(j,1)*x(1,i-j+1);
        end
    end
end

plot(sguji,'g');

%求出原始信号与估计信号的均方误差--------------------------------------
er=0;
for i=1:L
    er=er+((s(i)-sguji(i)).^2)/L;
end
disp(er);