% ARMA_GENE.M
%
% P. Flandrin, Sept. 25, 2003
%
% calcule et trace une r�alisation d'un mod�le ARMA
% d'ordres 0, 1, 2, 3 ou 4 + le spectre correspondant
%
% calls  : - ARMA_coeff.m

N = 1000; % nombre de points du processus
Nf = 1000; % nombre de points du spectre

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

n = 1; % num�ro de la figure de placement des p�les

f = linspace(0,.5,Nf);
z = exp(i*2*pi*f);

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

% pause

figure(n+1)

M = [ones(Nf,1) 1./z' 1./z'.^2 1./z'.^3 1./z'.^4]';

b = randn(1,N+5*p);
xx = filter(B,A,b);
x = xx(5*p+1:N+5*p);

subplot(311)
plot(x)
hold on
plot([1 N],[0 0])
title(['ARMA(',int2str(p),',',int2str(q),')'])
axis([1 N 1.1*min(x) 1.1*max(x)])
xlabel('temps')
set(gca,'YTick',[])
hold off

subplot(312)
S = (abs(B*M(1:q+1,:))).^2./(abs(A*M(1:p+1,:))).^2;
plot(f,log(S))
set(gca,'XTick',[])
ylabel('log(DSP)')

subplot(313)
spectrum(x)