% AR_IDENTIF_GC.M
%
% P. Flandrin, 30 sept. 2003 -- 14 nov. 2003
%
% exemple de modŽlisation AR du gain d'un mod�le ARMA
 
% ----- signal

N = 128; % nombre de points de la rŽponse impulsionnelle
Nfft = 512; % nombre de points pour le calcul de la FFT

p = 4; % ordre de la partie AR du mod�le (<= 4)
q = 0; % 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

c = randn(1,2*N);
Pth = mean((c-mean(c)).^2);
xx = filter(B,A,c);
x = xx(fix(N/2):fix(N/2)+N-1);

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

pause

xf = fft(x,Nfft)*sqrt(Nfft)/N;
subplot(222)
plot(log((abs(xf(1:Nfft/2))).^2),'r')
V = axis;
axis([1 Nfft/2 V(3) V(4)])
set(gca,'XTick',[])
xlabel('frequence')
ylabel('log(DSP)')

pause

% ----- identification

pmax = 15; % ordre max du mod�le estimŽ
cmax = 50; % nombre max de coeffs cepstraux

[AA,KK,CC,P,FPE,popt] = AR_levinson2(x,pmax,cmax);

figure(nf+1)
% clf

subplot(223)
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(224)
Nf = Nfft/2; % 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 = 1*(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',[])
W = axis;
axis([0 0.5 V(3) V(4)])
xlabel('frequence')
ylabel('log(DSP)')
title('pointilles = modele; trait plein = estimee') 
hold off

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

subplot(224)
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)]
% 
% [1 KK]
% 
% figure(nf+2)
% 
% %plot(0:cmax,[log(P(pp)) CC(pp,:)])
% 
% imagesc(abs(CC))
% colormap(jet)
% xlabel('indice cepstral')
% ylabel('ordre du mod�le')
% title('abs(coefficients cepstraux)')
% hold on
% plot([0 cmax+.5],[p-.5 p-.5],'w')
% plot([0 cmax+.5],[p+.5 p+.5],'w')
% hold off
% colorbar
