% AR_IDENTIF_RI.M
%
% P. Flandrin, 30 sept. 2003 -- 14 nov. 2003
%
% exemple de modˇlisation AR de la RI d'un mod¸le ARMA
 
% ----- signal

N = 50; % nombre de points de la rˇponse impulsionnelle

p = 4; % ordre de la partie AR du mod¸le (<= 4)
q = 2; % ordre de la partie MA du mod¸le (<= 4)

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

[A,B] = ARMA_coeff(p,q,nf);

pause

figure(nf+1)
clf

ind0 = 5;
imp = zeros(1,N);
imp(ind0) = 1;
x = filter(B,A,imp);

plot([1 N],[0 0],'k')
hold on
plot([ind0 ind0],[1.1*min(x) 1.1*max(x)],'k')
stem(x,'r')
title(['filtre ARMA(',int2str(p),',',int2str(q),')'])
axis([1 N 1.1*min(x) 1.1*max(x)])
xlabel('temps')
set(gca,'YTick',[])
set(gca,'XTick',[])
hold off

pause

% ----- identification

Pth = 1;

pmax = 20; % ordre max du mod¸le estimˇ

[AA,P,FPE,popt] = AR_levinson(x,pmax);

figure(nf+2)
clf

subplot(211)
plot(log(P),'b+-')
hold on
plot(log(FPE),'r+-')
plot(popt,log(FPE(popt)),'r','LineWidth',6)
xlabel('ordre du modele')
ylabel('log(erreur)')
title('bleu = puissance; rouge = FPE')
hold off

subplot(212)
Nf = 1000; % nombre de points du spectre
f = linspace(0,.5,Nf);
z = exp(i*2*pi*f);

M = [];
M(1,:) = ones(1,Nf);
for mm = 2:pmax+1
	
	M(mm,:) = 1./z.^(mm-1);
	
end

Aest = [1 AA(popt,1:popt)];
Sest = 1./(abs(Aest*M(1:popt+1,:))).^2;
S = Pth*(abs(B*M(1:q+1,:))).^2./(abs(A*M(1:p+1,:))).^2;
plot(f,log(S),'r--')
hold on
plot(f,log(Sest),'r','LineWidth',2)
set(gca,'XTick',[])
xlabel('frequence')
ylabel('log(DSP)')
title('pointilles = modele; trait plein = estimee') 
hold off

subplot(211)
hold on
[X,Y] = ginput(1);
pp = round(X(1));
plot(pp,log(P(pp)),'b','LineWidth',6)
hold off

subplot(212)
hold on
App = [1 AA(pp,1:pp)];
Spp = 1./(abs(App*M(1:pp+1,:))).^2;
plot(f,log(Spp),'b--','LineWidth',2)
hold off

A

[1 AA(pp,1:pp)]