% calcolo funzione densita' spettrale di potenza
% con il metodo del periodogramma mediato 
% per una realizzazione di p.a. AR(1) Gaussiano
% 
function [dsp,f] = periodmedar(N,M,a)

% in:  numero di campioni, N;
%      numero di segmenti, M (N deve essere divisibile per M,
%                             e la durata dei segm. N/M pari);
%      coefficiente di reazione, a;

% out: funzione densita' spettrale di potenza, dsp;
%      ascisse frequenze normalizzate, f;
%      uscita su video di grafico di periodogramma mediato.

sigmaw=sqrt(1-a^2);      % calcola la deviazione standard per il processo di forzamento
                         % per potenza unitaria del p.a. AR(1)
w=sigmaw*randn(1,N);     % genera una realiz. di durata N del processo di forzamento;
x=filter(1,[1 -a],w,0);  % filtraggio IIR del primo ordine della realizzazione
                         % (qui si trascura la presenza del transitorio).

lseg=N/M;                % lunghezza segmenti
x=reshape(x,lseg,M).';   % matrice di segmenti (uno per riga)

spettri=fft(x,[],2);              % calcola gli spettri complessi dei segmenti
dspsegm=abs(spettri).^2/(lseg-1); % calcola i periodogrammi dei segmenti
dsp=mean(dspsegm);                % calcola il periodogramma mediato

dsp=fftshift(dsp);             % centra la componente a frequenza zero
deltaf=1/lseg;                 % passo di campionamento in frequenza
f=[-lseg/2:lseg/2-1]*deltaf;   % scala delle frequenze

figure;                        % grafico dsp del processo
plot(f,dsp);
xlabel('f (norm.)');
ylabel('dsp');
v=axis;
axis([v(1) v(2) 0 v(4)])       % fissa a zero l'ordinata minima
