weinafilter_1.m

自适应维纳滤波器在matlab中的实现

clear all

%输入信号样本个数L，FIR滤波器阶数N%
L=input('L=');   
N=input('N=');
a=0.95;

%产生L个v(n)、u(n)、s(n)和x(n),在同一座标内绘出最后100个s(n)和x(n)
                         
v=1;               %信号方差
l=100;             %样本个数
u=0.1;            %步长
hi=-0.8;           %hi
h=zeros(l,1);
h(1)=0;            %h(0)
mp=0;              %平均超调量
ct=1;              %统计次数

for count=1:ct
    x=randn(1,l);
    meanx=sum(x)/l;
    covx=sum((x-meanx)*(x-meanx)')/l;
    x=0.5*(x-meanx)/sqrt(covx);
    w=randn(1,l);
    meanw=sum(w)/l;
    covw=sum((w-meanw)*(w-meanw)')/l;
    w=sqrt(v)*(w-meanw)/sqrt(covw);
    
     sn=hi*x;
    xn=s+w;

figure(1);                                %绘制s(n)(红色),x(n)(蓝色）
k=(L-99):1:L;                          
plot(k,s(k),'r',k,x(k),'k--');
legend('s(n)','x(n)',0);
title('s(n)和x(n)');
xlabel('n');ylabel('输入信号');


%利用L个v(n),s(n)，估计Rxx和rxs,调用矩阵求逆子程序计算h1，将N个理想的h(n)和估计的h1绘与同一座标内%
fxx=zeros(1,N);                        %计算fxx
for(i=1:N);
    for(k=0:(L-i-1));
        fxx(i)=(1/(L-i))*x(k+1)*x(k+i)+fxx(i);
    end;
end;
for(i=1:N);                            %生成Rxx矩阵
    for(k=1:N);
        Rxx(k,i)=fxx(abs(k-i)+1);
    end;
end;

rxs=zeros(1,N);                        %计算rxs
for(i=1:N);
    for(k=0:(L-i-1));
        rxs(i)=(1/(L-i))*x(k+1)*s(k+i)+rxs(i);
    end;
end;

h1=(inv(Rxx))*rxs';                   

for(i=1:N);                           %绘制h(理想滤波器)，h1(估计滤波器)
    h(i)=0.238*(0.724)^i*u(i);
end;

figure;
k=1:1:N;
plot(k,h(k),'r',k,h1(k),'k--');
title('理想的h(n)和估算的h(n)');
legend('理想的h(n)','估算的h(n)');
xlabel('n');

%重新生成L点随机序列%
w=-sqrt(3*(1-(0.95^2)))+2*sqrt(3*(1-(0.95^2)))*rand(1,L);                        
v=-sqrt(3)+2*sqrt(3)*rand(1,L);
u=ones(1,L);

s(1)=1;                                 %计算s(n)
i=2;
while(i<=L);
   s(i)=a*s(i-1)+w(i);
   i=i+1;
end;

i=1;                                   %计算x(n)
while(i<=L);
    x(i)=s(i)+v(i);
    i=i+1;
end;

%利用公式(1－15)，进行理想维纳滤波得L个SI(n),将最后100个S(n)和SI(n)绘制于同一坐标内%
S=conv(h,x);
SI(1)=S(1);

for(i=2:L);
    SI(i)=0.724*SI(i-1)+0.238*x(i);
end;

figure;
k=(L-99):1:L;
plot(k,s(k),'r',k,SI(k),'k--');
title('s(n)和Si(n)');
legend('s(n)','Si(n)',0);
xlabel('n');ylabel('理想滤波输出');

%由h1(n),根据公式(1-6)对x(n)进行过滤得L个SR(n)将最后100个s(n)和SR(n)绘制与同一坐标内%
SR=conv(h1,x);
figure;
k=(L-99):1:L;
plot(k,s(k),'r',k,SR(k),'k--');
title('s(n)和Sr(n)');
legend('s(n)','Sr(n)',0);
xlabel('n');ylabel('实际滤波输出');

%由L个x(n),s(n),SI(n),SR(n),根据公式(1-17)，统计并打印出EX2,EI2,ER2%
EX2=0;
EI2=0;
ER2=0;

for(i=1:L);
    EX2=1/L*((x(i)-s(i))^2)+EX2;               %计算接受信号与原信号的误差
end;

for(i=1:L);
    EI2=1/L*((SI(i)-s(i))^2)+EI2;               %计算维纳滤波器理想滤波
end;

for(i=1:L);                                  %计算N点的维纳滤波器理想滤波
    EI21=1/L*((S(i)-s(i))^2)+EI2;
end;

for(i=1:L);                                   %计算估计滤波器滤波
    ER2=1/L*((SR(i)-s(i))^2)+ER2;
end;

EX2
EI2
ER2