% 本实验中利用近似方法，即最佳FIR维纳滤波方法，在计算机上实现随机信号的维纳滤波。

% w(n)是零均值,方差为(1-a^2)的均匀分布白噪声

% s(n)为真实信号：s(n)=a*s(n-1)+w(n)

% v(n)是与s(n)互不相关的均匀分布白噪声，其均值为零，方差为1

% x(n)为接收到的添加了白噪声的信号：x = s + v

% si(n) 为使用理想维纳滤波器滤波得到结果

% y(n) 为使用近似FIR维纳滤波器滤波得到的结果

% If you have any problem, please email me:

qijiaxing@live.cn

% This program is written by QiJiaxing

% 2009.11.20

clc

clear all

L = 500;  % signal length

N = 20;  % length of the FIR

filter

a = 0.95;

% white noise with mean of 0 and var of (1-a^2)

w =  sqrt( 12 * (1 - a^2)) * ( rand(1,L) -

0.5 );

% true signal: s(n)=a*s(n-1)+w(n)

s = zeros(1,L); s(1) = w(1);

for ii = 2:L

s(ii) =

a * s(ii-1) + w(ii);

end

% white noise with mean of 0 and var of 1

v =  sqrt( 12 ) * ( rand(1,L) - 0.5 );

% received signal: x = s + v

x = s + v;

%  r_xx is the autocorrelation of x

r_xx = xcorr( x );

% R_xx is the N-dimentional autocorrelation matrix of x

R_xx = zeros( N );

for ii = 1 : N

R_xx( :

, ii ) = r_xx( L+1-ii : L+1-ii+N-1 )';

end

%  r_xs is the cross-correlation of x and

s

r_xs = xcorr( x , s );

r_xs = r_xs( L : L+N-1 )';

% according to R_xx * h_FIR = R_xs

% we can calculate the h_FIR

h_FIR = inv(R_xx) * r_xs;

h_FIR = h_FIR';

% y is the reslut signal filtered by h_FIR

y = cconv( h_FIR , x, L );

% si is the result signal filtered by ideal Wiener filter

h

n = 0:L-1;

h = 0.238 * ( 0.724 .^ n );

si  = cconv( h , x, L );

%% plot the true signal s(n) and the received signal x(n) of

the last 100

% points

t = L - 99 : L;

figure(1);

plot( t , x(t), '--b' , t , s(t) , 'r' );

legend( 'x(n)' , 's(n)' , 0 );

title('The true signal s(n) and the received signal

x(n)');

xlabel('n');ylabel('Signal Amplitude');

%% plot the ideal impulze response h(n) and the FIR impulze

response

% h_FIR(n) of N points

n = 1:N;

figure(2);

plot( n,h(n),'--b',n,h_FIR(n),'r' );

legend('The Ideal Wiener Filter h(n)','Approximate FIR Filter

h\_FIR(n)',0);

title('The Impulze Response of Ideal Wiener Filter h(n) and

the Impulze Response h\_FIR(n)');

xlabel('N points');ylabel('Amplitude');

%% plot the true signal s(n) and the result signal y(n)

filtered by

% h_FIR

figure(3);

plot( t, s(t), '--b', t, y(t), 'r' );

legend('True Signal s(n)','Result Signal y(n)',0);

title('The True Signal s(n) and The Result Signal y(n) Filter

by h\_FIR');

xlabel('n');ylabel('Signal Amplitude');

%% plot the true signal s(n) and the result signal si(n)

filtered by h

figure(4);

plot( t, s(t), '--b', t, si(t), 'r' );

legend('True Signal s(n)','Result Signal si(n)',0);

title('The True Signal s(n) and The Result Signal si(n)

Filtered by h');

xlabel('n');ylabel('Signal Amplitude');