% Analisi spettrale di segnali stazionari: due sinusoidi con rumore
% Uso della funzione "periodogram" per la stima dello spettro di potenza
% Uso della funzione pwelch per la stima della PSD attraverso lo "spettrogramma mediato"
clear all;

%  somma di due sinusoidi con rumore
Fs= 1000;	% [Hz] frequenza di campionamento
L= 512;	% durata segnale (# campioni)
t= (0:L-1)./Fs;	% ascisse campioni
A0= 1;  A1= 0.75;	% Ampiezze sinusoidi
% i valori seguenti corrispondono al caso visto a lezione che creava problemi
f0= Fs/14;  f1= Fs/7.5;	% frequenze non normalizzate sinusoidi
% i valori seguenti corrispondono al caso  fortuito visto a lezione che individua
% tutti campioni nulli della DFT (v periodico con periodo L)
% f0= Fs/16;  f1= Fs/8;	% pulsazioni normalizzate sinusoidi
% rumore Gaussiano a componenti incorrelate con potenza media sigma^2
randn('state',0);
sigma= 1; wt= sigma*randn(1,L);
% segnale campionato
vt= A0*cos(2*pi*f0*t)+A1*cos(2*pi*f1*t)+wt;
% tempo
figure(1); clf;
stem(t,vt,'r');
xlabel('t [s]'); ylabel('v(t)');
axis([0 (L-1)/Fs -4 4]);
title('due sinusoidi + rumore');
% frequenza
fvec= linspace(-Fs/2,Fs/2,length(t));
Vf= abs(fft(vt));
figure(2); clf;
plot(fvec,fftshift(Vf),'r');
xlabel('f [Hz]'); ylabel('|V(f)|');
axis([-Fs/2 Fs/2 0 L/2]);
title('spettro (DFT) del segnale+rumore');

% Valutiamo ora la densità spettrale di potenza del segnale (aleatorio) vt
% Periodogramma: deifinito come |V(f)|^2/(L*Fs)
figure(5); clf;
nfft= 1024;
[Pvv,fv]= periodogram(vt,[],'twosided',nfft,Fs);
% dove l'argomento è vuoto [], si può usare il nome di una finestra
% NOTA: in matlab, la finestra viene normalizzata in ampiezza in modo da non alterare il calcolo della pot. media
% [Pvv,fv]= periodogram(vt,hamming(L),'twosided',nfft,Fs);
plot(fv,Pvv,'m');
xlabel('f [Hz]'); ylabel('P_V(f)');
title(' Densità spettrale di potenza (periodogramma)');
% Calcolo della potenza media
sprintf('Potenza media segnale+rumore: teorico= %f , stimato= %f',A0^2/2+A1^2/2+sigma^2,(Fs/length(Pvv))*sum(Pvv))
% NOTA: Fs/n, dove n è la dimensione do Pvv, è il "df" con cui integro la densità Pvv (approssimata con rettangoloide)

% si può visualizzare direttamente in scala logaritmica col comando
figure(6); clf;
periodogram(vt,[],'twosided',nfft,Fs);

% NOTA IMPORTANTE: anche aumentando la lunghezza del segnale, le oscillazioni (statistiche) dello spettro del rumore non si appiattiscono

% Un metodo per aumentare il rapporto segnale-rumore nello spettro è quello di usare il "metodo di Welch"
Lfin= 32;	% default: 8 finestre con 50% overlap
finestra= hamming(Lfin);	% default =hamming()
Loverlap= 8;	% default =Lfin/2
nfftWelch= 1024;	% default = 256 (o 512, 1024...) > lunghezza segnale
figure(7); clf;
pwelch(vt,finestra,Loverlap,nfftWelch,Fs,'twosided');

% NOTA: lo spettro del rumore viene "ammorbidito" e attenuato dalla media
%    ...e questo a spese della Risoluzione, diminuita sia da finestra corta (Lfin) che da finestre sagomate (es. hamming)

% chiamata così, la stima di Welch coincide col periodogramma semplice
figure(8); clf;
pwelch(vt,length(vt),0,nfft,Fs,'twosided');