filter.m

数字信号处理中

clear all
close all

fprintf('\n\n下面开始实验： \n');
L=input('请输入信号的样本个数，L=');
N=input('请输入FIR滤波器的阶数，N=');

a=0.95;
dw=1-a*a;

w = -sqrt(3*dw) + 2*sqrt(3*dw) * rand(1,L);
v = -sqrt(3)+2*sqrt(3)*rand(1,L);

s=zeros(1,L);
x=zeros(1,L);
for i=2:L
    s(i)=a*s(i-1)+w(i); 
end
x=s+v;

plot(x(1,[(L-99):L]),'r');
grid on;
hold on;
title('x(n)和s(n)比较');
plot(s(1,[(L-99):L]),'--b');
legend('x(n)','s(n)');
xlabel('采样点');
ylabel('幅值');

%************************************************************************
% 求CORxx(n)和CORxs(n)
CORxx=zeros(1,N);
CORxs=zeros(1,N);
for m=1:N,
    for i=1:L-(m-1),
        CORxx(m)=CORxx(m)+x(i)*x(i+(m-1));
        CORxs(m)=CORxs(m)+x(i)*s(i+(m-1));
    end
    CORxx(m)=CORxx(m)/(L-(m-1));
    CORxs(m)=CORxs(m)/(L-(m-1));
end

% 给Rxx(n)赋值
Rxx=zeros(N,N);
for i=1:N,
    for j=1:N,
        t=abs(i-j);
        Rxx(i,j)=CORxx(t+1);
    end
end

% 给rxs(n)赋值
rxs=zeros(N,1);
for i=1:N;
    rxs(i,1)=CORxs(i);
end

% 求估计的h(n)
h_opt=(Rxx^(-1)*rxs);

% 理想的h(n)
u=ones(N,1);
for n=1:N
    h(n)=0.238*(0.724)^n*u(n);
end 

% 画图
a=figure;
plot(h_opt,'r');
grid on;
hold on;
title('估计值h_opt（n)和理想值h(n)比较');
plot(h,'--b');
legend('h_opt(n)','h(n)');

%*************************************************************************
% 求解理想维纳滤波的结果sI(n)
sI=zeros(1,L);
for i=2:L
    sI(i)=0.724*sI(i-1)+0.238*x(i);
end
b=figure;
plot(sI(1,[(L-99):L]),'r');
grid on;
hold on;
title('信号s(n)和理论sI(n)比较');
plot(s(1,[(L-99):L]),'--b');
legend('sI','s');

%*************************************************************************
% 求解实际维纳滤波的结果sR(n)
sR = conv(x,h_opt);
c=figure;
plot(sR(1,[(L-99):L]),'r');
grid on;
hold on;
title('信号s(n)和实际滤波值sR(n)比较');
plot(s(1,[(L-99):L]),'--b');
legend('sR','s');

%*************************************************************************
% 计算误差
E_sI=0;
E_sR=0;
E_x=0;
for i=1:L,
    E_x=E_x+(x(i)-s(i))^2;
    E_sI=E_sI+(sI(i)-s(i))^2;
    E_sR=E_sR+(sR(i)-s(i))^2;
end
ex=E_x/L;
eI=E_sI/L;
eR=E_sR/L;
fprintf('不进行滤波的误差为：%f \n',ex);
fprintf('维纳滤波理论误差为：%f \n',eI);
fprintf('FIR逼近维纳滤波的实际误差为：%f \n',eR);
%**************************************************************************