wienerfilter.m

这是一个wiener滤波程序

function y=wienerfilter(x,Rxx,Rxd,N)
%进行维纳滤波
%x是输入信号，Rxx是输入信号的自相关向量
%Rxx是输入信号和理想信号的的互相关向量，N是维纳滤波器的长度
%输出y是输入信号通过维纳滤波器进行维纳滤波后的输出
h=yulewalker(Rxx,Rxd,N);%求解维纳滤波器系数
t=conv(x,h);            %进行滤波
Lh=length(h);           %得到滤波器的长度
Lx=length(x);           %得到输入信号的长度
y=t(double(uint16(Lh/2)):Lx+double(uint16(Lh/2))-1);%输出序列y的长度和输入序列x的长度相同


%以下是维纳滤波器系数的求解
function h=yulewalker(A,B,M)   
%求解Yule-Walker方程
%A是接收信号的自相关向量为 Rxx(0),Rxx(1),......,Rxx(M-1)
%B是接收信号和没有噪声干扰信号的互相关向量为 Rxd(0),Rxd(1),......,Rxd(M-1)
%M是滤波器的长度
%h保存滤波器的系数
%例如A=[6,5,4,3,2,1];B=[100,90,120,50,80,200];M=6;
%求解出h=[26.4286  -20.0000   50.0000  -50.0000  -45.0000   81.4286]'
T1=zeros(1,M);%T1存放中间方程的解向量
T2=zeros(1,M);%T2存放中间方程的解向量
T1(1)=B(1)/A(1);
T2(1)=A(2)/A(1);
X=zeros(1,M);
for i=2:M-1
temp1=0;
temp2=0;
    for j=1:i-1
        temp1=temp1+A(i-j+1)*T1(j);
        temp2=temp2+A(i-j+1)*T2(j);
    end
    X(i)=(B(i)-temp1)/(A(1)-temp2);
    for j=1:i-1
        X(j)=T1(j)-X(i)*T2(j);
    end
    for j=1:i
        T1(j)=X(j);
    end
temp1=0;
temp2=0;
    for j=1:i-1
        temp1=temp1+A(j+1)*T2(j);
        temp2=temp2+A(j+1)*T2(i-j);
    end
    X(1)=(A(i+1)-temp1)/(A(1)-temp2);
    for j=2:i
        X(j)=T2(j-1)-X(1)*T2(i-j+1);
    end
    for j=1:i
        T2(j)=X(j);
    end
end
temp1=0;
temp2=0;
for j=1:M-1
	temp1=temp1+A(M-j+1)*T1(j);
	temp2=temp2+A(M-j+1)*T2(j);
end
X(M)=(B(M)-temp1)/(A(1)-temp2);
for j=1:M-1
	X(j)=T1(j)-X(M)*T2(j);
end
h=X;