% ARMA_WIENER.M
%
% P. Flandrin, Sept. 25, 2003
%
% calcule et trace la filtrˇe de Wiener d'une
% rˇalisation d'un mod¸le ARMA (d'ordres 0, 1, 2, 3 ou 4)
% noyˇ dans un bruit ARMA
%
% calls  : - ARMA_coeff.m

N = 200; % nombre de points du processus
Nout = 100;
Nf = 1024; % nombre de points du spectre
f = linspace(0,.5,Nf);
z = exp(i*2*pi*f);
M = [ones(Nf,1) 1./z' 1./z'.^2 1./z'.^3 1./z'.^4]';

ps = 4; % ordre de la partie AR du mod¸le de signal (<= 4)
qs = 0; % ordre de la partie MA du mod¸le de signal (<= 4)

ns = 1; % numˇro de la figure de placement des p™les

pb = 2; % ordre de la partie AR du mod¸le de bruit (<= 4)
qb = 2; % ordre de la partie MA du mod¸le de bruit (<= 4)

nb = ns+1; % numˇro de la figure de placement des p™les

SNR = 0;
alpha = 10^(-SNR/20);

% ----- signal

[As,Bs] = ARMA_coeff(ps,qs,ns);
bs = randn(1,N);
s = filter(Bs,As,bs);
s = s/sqrt(sum(s.^2));

% pause

% ----- bruit

[Ab,Bb] = ARMA_coeff(pb,qb,nb);
bb = randn(1,N);
b = filter(Bb,Ab,bb);
b = b/sqrt(sum(b.^2));

y = s + alpha*b;

% ----- filtre

Ss = (abs(Bs*M(1:qs+1,:))).^2./(abs(As*M(1:ps+1,:))).^2;
Sb = (abs(Bb*M(1:qb+1,:))).^2./(abs(Ab*M(1:pb+1,:))).^2*alpha^2;

Hf = Ss./(Ss+Sb);

figure(nb+1)
clf

subplot(211)

plot(f,10*log10(Ss),'r-')
hold on
plot(f,10*log10(Sb),'b--')
plot(f,10*log10(Ss+Sb),'b','LineWidth',2)
plot(f,10*log10(Hf),'k','LineWidth',2)
xlabel('frˇquence')
ylabel('DSP (dB)')
title('bleu continu = observation; bleu tirets = bruit; rouge = signal; noir = filtre de Wiener')
hold off

pause

H = zeros(1,2*Nf);
H(1:Nf) = Hf;
H(Nf+2:2*Nf) = Hf(Nf:-1:2);

h = real(ifft(H));

ind = find(abs(h(1:Nf))>=h(1)/1000);
B = h(ind);
LB = length(B);
BB = [fliplr(B(2:LB)) B];
dec = LB-1;

sest = filter(BB,1,y);

SNRout = 10*log10(sum(s(1:Nout).^2)/sum((s(1:Nout)-sest(1+dec:Nout+dec)).^2));

subplot(212)

plot(s(1:Nout),'r+-')
hold on
plot(y(1:Nout),'b+-')
plot(sest(1+dec:Nout+dec),'k+-')
xlabel('temps')
title(['bleu continu = observation; rouge = signal original; noir = signal estime; SNRin = ',int2str(SNR),' ; SNRout = ',int2str(SNRout)])
hold off
zoom on