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);
%**************************************************************************