dsp1.m

使用维纳滤波器对信号进行处理

L=input('请输入信号样本个数L:');
N=input('请输入FIR滤波器阶数N:');
Ex=0;                               %x(n)对s(i)的均方误差
Ei=0;                               %si(n)对s(i)的均方误差
Er=0;                               %sr(n)对s(i)的均方误差
b1=sqrt(3);                         %产生v(n)的参数
b2=0.5408;                          %产生w(n)的参数
a=0.95;                             %设定a值

v=2*b1*rand(1,L)-b1;                %产生v(n)，方差为1的随机数
u = zeros(1,L);
for i = 1:L,                        %产生u(n)
    u(i) = 1;
end
w=2*b2*rand(1,L)-b2;                %产生w(n)，方差为1－a^2
s = zeros(1,L);
s(1)=0;                             %初始化s(1)=0
for i = 2:L,                        %得到信号s(n)
    s(i) = a*s(i-1)+w(i);
end
x = zeros(1,L);
x = s + v;                          %得到含有噪音的信号x(n)
n=L-100:L;

plot(n,x(n),'k:',n,s(n),'b');       %在同一坐标系中画出最后100个s(n),%x(n)
legend('x(n)','s(n)');              %标注波形的归属
title('信号s(n)和噪声信号混合x(n)'); %加图形标题
xlabel('n');                        %加X轴说明
ylabel('x(n)--s(n)');               %加Y轴说明

Rxx = zeros(N,N);                   %Rxx为x(n)的自相关
for i = 1:N,                        %i为行数
    for j = 1:N,                    %j为列数
        teR=0;   
        te=abs(i-j);                %括号中的值是行数减去列数
        for k = 1:L-te,             %由式1-18得到Rxx 
            teR=teR+x(k)*x(k+te);
        end
        teR=teR/(L-te);
        Rxx(i,j)=teR;               %赋值给Rxx
    end
end
rxs = zeros(N,1);                   %rxs为x(n)和s(n)的相关数
for m = 0:N-1,                      %由式1-19得到
    ter=0;
    for i = 1:L-m,
        ter=ter+x(i)*s(i+m);
    end
    ter=ter/(L-m);
    rxs(m+1)=ter;                     
end
h_e=Rxx^(-1)*rxs;                    %由式1-11求出h_e
h=zeros(N,1);                        %维纳滤波器理想脉冲相应
for i = 1:N,
    h(i)=0.238*0.724^i*u(i);         %求出理想h
end
n = 1:N;

figure;
plot(n,h_e(n),':*',n,h(n),'-o');       %在同一坐标系中绘出h和它的估计值
legend('h_e(n)','h(n)');
title('估计h_e(n)和理想h(n)');        %加图形标题
xlabel('n');                         %加X轴说明
ylabel('h_e(n)--h(n)');              %加Y轴说明
v=2*b1*rand(1,L)-b1;                 %产生v(n)，方差为1的随机数
w=2*b2*rand(1,L)-b2;                %产生w(n)，方差为1－a^2
s = zeros(1,L);
s(1)=0;                             %初始化s(1)=0
for i = 2:L,                          %得到信号s(n)
    s(i) = a*s(i-1)+w(i);
end
x = zeros(1,L);
x = s + v;                           %得到含有噪音的信号x(n)

si_e=zeros(1,L);
si_e(1)=0;                           
for i = 2:L,                         %由1-15得到si_e(n)
    si_e(i)=0.724*si_e(i-1)+0.238*x(i);
end
n = L-100:L;

figure;
plot(n,si_e(n),':',n,s(n));          %将s与理想维纳滤波的sl值绘于同一坐标系中
legend('si_e(n)','s(n)');
title('si_e(n)和s(n)');              %加图形标题
xlabel('n');                         %加X轴说明
ylabel('si_e(n)--s(n)');             %加Y轴说明

sr_e=zeros(1,L);
for i = N:L,
    for m = 0:N-1,
        sr_e(i)=sr_e(i)+h_e(m+1)*x(i-m);
    end
end                                   %由(1-6)算出est_sr(n)
n = L-100:L;

figure;
plot(n,sr_e(n),':',n,s(n));           %将s与S_R值绘于同一坐标系中
legend('sr_e(n)','s(n)');
title('sr_e(n)和s(n)');               %加图形标题
xlabel('n');                          %加X轴说明
ylabel('sr_e(n)--s(n)');              %加Y轴说明

for i = 1:L,
    Ex=Ex+(x(i)-s(i))^2;
end
Ex=Ex/L,
for i = 1:L,
    Ei=Ei+(si_e(i)-s(i))^2;
end
Ei=Ei/L,
for i = 1:L,
    Er=Er+(sr_e(i)-s(i))^2;
end
Er=Er/L,
figure;
plot(1,Ex,'p',1,Ei,'*',1,Er,'o');     %将EX^2,EL^2,ER^2,绘于同一坐标系中
legend('Ex','Ei','Er');
title('Ex,Ei,Er');                    %加图形标题