Tc = 0.01;
nper = 4;  % numero di periodi considerati in ciascuna sotto-sequenza
N = 1/Tc * nper;
nsseq = 50; % numero di sotto-sequenze 
t = 0:Tc:nper*nsseq-Tc;
si = sin(2*pi*t) + sin(2*pi*3*t) + sin(2*pi*5*t) + 6*randn(1,length(t));
figure(1); plot(t,si,'b-'); axis([0 (nper*nsseq)/10 -8 8]), title('segnale con rumore'); pause,
%
[I,f] = psd_Welch(si,Tc,nsseq);
figure(2), plot(f, I(1:N/2), 'b-'), 
           xlabel('frequenza [Hz]'),
           title(['numero di sotto-sequenze = ' num2str(nsseq)]),
           pause,
           axis([0 10 20 140]); grid minor;
%
% NB: confrontare la bontà di questo risultato, che considera la sequenza come segnale casuale e
% dunque ne valuta lo spettro di potenza, rispetto a quello ottenuto alla Lezione 3, in cui si considerava 
% la sequenza come segnale deterministico e si calcolava semplicemente lo spettro. 
%

%
% Verifica sperimentale della relazione di Parseval, e cioè che l'integrale (sommatoria) dello 
% spettro di potenza è una misura della potenza media del segnale:
%
disp(['valore quadratico medio del segnale = ' num2str(sum(sum(si.^2))/(N*nsseq))]);
disp(['valor medio della stima dello spettro di potenza = ' num2str(sum(I)/N)]);



%%%
%%%  Calcoliamo in vari modi la risposta in frequenza del sistema meccanico con le tre masse:
%%%
N = 1024*10;
%
% costruisco il modello continuo del sistema meccanico:
[A,B,C,D]=build_sistema_meccanico_3gdl();
%
% definisco il periodo di campionamento:
Ts = 0.0001;
%
%
% calcolo analitico della risposta in frequenza:
%
if 1  % possiamo eseguire questa sezione di codice solo se abbiamo il "control system toolbox"
  ss_d_3gdl = c2d(ss(A,B(:,1),C,D(:,1)), Ts);
  tf_d_3gdl = tf(ss_d_3gdl);
  [b,a] = tfdata(tf_d_3gdl,'v');
  figure(1); clf, freqz(b,a);
end
%
%
%  Stima della risposta in frequenza da una sequenza di dati "sperimentali":
%
% 1. risposta ad un ingresso sinusoidale:
omega_nu = (2*pi/N) * 0.1*(N/2);
u = cos(omega_nu * (1:N)) + 0.1*randn(1,N);
y = simula_sistema_meccanico_3gdl(A,B,C,D,u,Ts);
figure(2); clf, hold on, plot(1:N/4,u(1:N/4),'g-', 1:N/4,y(1:N/4)*1e7,'b-'); title('"y" è amplificato !'); pause,
[Puu, f] = psd(u);
[Pyy, f] = psd(y);
figure(1); subplot(2,1,1); hold on, plot(f, 20*log10( sqrt(Pyy / Puu) ), 'r--'); 
           title('stima (non parametrica) della risposta in frequenza'); 
subplot(2,1,1); axis([0 0.2 -250 -140]);
subplot(2,1,2); axis([0 0.2 -250 0]);
pause
% --> notare che la stima ha validità nello spettro di frequenze del segnale di ingresso, che qui è una sinusoide pura
%
% 3. risposta ad un ingresso aleatorio a largo spettro, es. un rumore bianco:
u = 1 * randn(N,1);
y = simula_sistema_meccanico_3gdl(A,B,C,D,u,Ts);
figure(3); clf, hold on, plot(1:N/4,u(1:N/4),'g-', 1:N/4,y(1:N/4)*1e7,'b-'); title('"y" è amplificato !'); pause,
[Puu, f] = psd(u);
[Pyy, f] = psd(y);
figure(1); subplot(2,1,1); hold on, plot(f, 20*log10( sqrt(Pyy / Puu) ), 'm--'); 
           title('stima (non parametrica) della risposta in frequenza'); 

pause



%
% test di Bartlett (o del periodogramma cumulato): serve per verificare se una sequenza è rumore bianco:
%
N = 4096*5;
var_s = 6;
s = randn(N,1) * sqrt(var_s);
figure(3); plot(s); title('rumore bianco'); pause;
[CPss_wn,f] = testBartlett_wn(s);
figure(4); plot(f, CPss_wn); title('periodogramma cumulato per un rumore bianco'); print -depsc LEZ_4_figure_4.epsc; pause,
s = randn(N,1) * sqrt(var_s) + sin(2*pi/(N/400)*(1:N))';
figure(5); plot(s); title('rumore bianco più componente sinusoidale'); pause;
[CPss,f] = testBartlett_wn(s);
figure(5); clf, plot(f, CPss_wn, 'r--'); hold on,
           plot(f, CPss); title('periodogramma cumulato per un rumore bianco più componente sinusoidale');
pause

