%
% utilizzo in pratica della DFT
%

echo on,
%
% troncamento di una sequenza sinusoidale pura ad un istante che non è un multiplo del periodo della
% sinusoide:
%
fc = 1;
Tc = 1/fc;
N = 256;
t = 0:Tc:(N-1)*Tc;
f = 0.0156;
s = sin(2*pi*f*t);   % segnale sinusoidale alla frequenza "f" (Hz)
% FFT di un multiplo intero del periodo:
N = 256;   % NB: round(fc/f) è circa = 64 , quindi si considerano 4 periodi esatti !!!   
t = 0:Tc:(N-1)*Tc;
s = sin(2*pi*f*t);
S = fft(s);
figure(3); clf,
subplot(2,1,1); plot(t, s, 'b-'); title('sequenza lunga un multiplo intero del periodo');
subplot(2,1,2); plot(0:1/(Tc*N):1/Tc-1/N, abs(S), 'b-o'); axis([0 0.3 0 max(abs(S))*1.1]);
pause
% FFT di un multiplo NON intero del periodo:
%N = 240; 
%t = 0:Tc:(N-1)*Tc;
tail = [zeros(1,240) s(N-[15:-1:0])];
sc = s - tail;
Sc = fft(sc);
figure(30); clf,
subplot(2,1,1); plot(t, sc, 'b-o'); title('sequenza lunga un multiplo NON-intero del periodo'); hold on,
subplot(2,1,2); plot(0:1/(Tc*N):1/Tc-1/N, abs(Sc), 'b-o'); axis([0 0.3 0 max(abs(S))]);  hold on,
pause
subplot(2,1,1); plot(t, -tail, 'r-+');
pause
subplot(2,1,2); plot(0:1/(Tc*N):1/Tc-1/N, abs(fft(tail)), 'r-*');
                
% notare che il prolungamento per periodicità introduce un "salto", discontinuità il cui spettro è a
% larga banda; questo spiega il risultato.
pause
% rimedio: FFT di un multiplo NON intero del periodo ma CON FINESTRATURA OPPORTUNA DEI DATI:
N = 240;
t = 0:Tc:(N-1)*Tc;
w = hamming(N); figure(31); plot(w,'b-o'); title('finestra di Hamming'); %print -depsc LEZ_3_figure_31.epsc; pause,
s = sin(2*pi*f*t) .* w;
S = fft(s);
figure(32); clf, 
subplot(2,1,1); plot(t, s, 'b-'); title('sequenza con finestratura di Hamming');
subplot(2,1,2); plot(0:1/(Tc*N):1/Tc-1/N, abs(S), 'b-o'); axis([0 0.3 0 max(abs(S))]);
%print -depsc LEZ_3_figure_32.epsc; 
pause
%
% DFT della finestra rettangolare e di quella di Hamming
%
N = 256;
w = [hamming(7*N/8) zeros(1,N/8)];
r = [ones(1,7*N/8) zeros(1,N/8)];  % finestra rettangolare
WdB = 20*log10( abs( fft(w) ) );
RdB = 20*log10( abs( fft(r) ) );
figure(33); clf, hold on,
            semilogy(0:1/(Tc*N):1/Tc-1/N, WdB, 'b-');
            semilogy(0:1/(Tc*N):1/Tc-1/N, RdB, 'r-'); 
            axis([0 0.06 -50 50]);
            legend('Hamming','rettangolare')
pause


%
% effetto "picket-fence" e tecnica di "zero-padding" (rimedio):
%
Tc = 0.1;  
N = 256;
t = 0:Tc:(N-1)*Tc;
f = 1/(Tc*N) * 12   % la frequenza cade ad una frequenza discreta risoluta dalla DFT
s = sin(2*pi*f*t);
S = fft(s);
figure(4); clf, 
subplot(2,1,1); plot(t, s, 'b-'); title('sinusoide a frequenza coincidente con una frequenza discreta della TFD');
subplot(2,1,2); plot(0:1/(Tc*N):1/Tc-1/N, abs(S), 'b-o'); axis([0 2 0 max(abs(S))]);
                xlabel('frequency [Hz]'); title(['risoluzione in frequenza = ' num2str(1/(Tc*N)) ' [Hz]']);
pause,
f = 1/(Tc*N) * 12.5   % la frequenza cade giusto a metà tra due frequenze discrete della DFT
s = sin(2*pi*f*t);
S = fft(s);
figure(40); clf, 
subplot(2,1,1); plot(t, s, 'b-'); title('sinusoide a metà tra due frequenze discrete della TFD');
subplot(2,1,2); plot(0:1/(Tc*N):1/Tc-1/N, abs(S), 'b-o'); axis([0 2 0 max(abs(S))]);
                xlabel('frequency [Hz]'); title(['risoluzione in frequenza = ' num2str(1/(Tc*N)) ' [Hz]']);
                % NB: ci sono due valori consecutivi dello spettro che indicano la sinusoide;
pause,
% rimedio: "zero-padding" (cioè aggiunta di zeri alla sequenza, previa finestratura per non creare
% discontinuità
f = 1/(Tc*N) * 12.5   % la frequenza cade giusto a metà tra due frequenze discrete della DFT
w = hamming(N); 
s = sin(2*pi*f*t) .* w;
s = [s zeros(1,N)];
N = length(s);
t = 0:Tc:(N-1)*Tc;
S = fft(s);
figure(40); clf, 
subplot(2,1,1); plot(t, s, 'b-'); title('sinusoide a metà tra due frequenze discrete della TFD, con ZERO-PADDING');
subplot(2,1,2); plot(0:1/(Tc*N):1/Tc-1/N, abs(S), 'b-o'); axis([0 2 0 max(abs(S))]);
                xlabel('frequency [Hz]'); title(['risoluzione in frequenza = ' num2str(1/(Tc*N)) ' [Hz]']);
pause,



%
% risoluzione in frequenza: aumentando il numero di campioni, a parità di frequenza di campionamento, aumenta la risoluzione in frequenza
%
Tc = 0.1;  
N1 = 256;
N = N1;
t = 0:Tc:(N-1)*Tc;
s1 = sin(2*pi*1*t) + sin(2*pi*1.02*t) + sin(2*pi*1.07*t);
for i=1:5
  N = N1*i;
  t = 0:Tc:(N-1)*Tc;
  s = sin(2*pi*1*t) + sin(2*pi*1.02*t) + sin(2*pi*1.07*t);
  S = fft(s);
  S1 = fft(s1);
  figure(5), clf, hold on,
             plot(0:1/(Tc*N):1/Tc-1/N, abs(S), 'b-'), 
             plot(0:1/(Tc*N):1/Tc-1/N, abs(S1), 'm--'), 
             axis([0 2 0 max(abs(S))]);
             xlabel('frequency [Hz]'); title(['risoluzione in frequenza = ' num2str(1/(Tc*N)) ' [Hz]']); pause,  
  s1 = [s1 zeros(1,N1)];             
end;




%
% presenza di rumore bianco a media nulla nei dati:
%
Tc = 0.01;
N = 1/Tc;
nper = 200;
t = 0:Tc:nper-Tc;
s = sin(2*pi*t) + sin(2*pi*3*t) + sin(2*pi*5*t);
echo off
figure(1); plot(t,s,'b-'); axis([0 nper/10 -8 8]), title('segnale senza rumore'); pause,
echo on,
% aggiungiamo il rumore:
s = s + 6*randn(1,length(s));
echo off
figure(1); plot(t,s,'b-'); axis([0 nper/10 -8 8]), title('segnale con rumore'); pause,
echo on
% vediamo la FFT mediata su più periodi:
echo off
s = reshape(s,N,nper);
for i=1:10:nper
  figure(2), plot(0:N-1, mean(abs(fft(s(:,1:i))),2), 'b-'), title(['numero di periodi = ' num2str(i)]), pause,
end



