%
% comportamento della DFT per segnali non-stazionari:
%
N = 512;
omega_nu = (2*pi/N) * 0.3*(N/2);
% segnale stazionario: 
u = (sin(omega_nu * (1:N)) + sin(1.5*omega_nu * (1:N)) + sin(2*omega_nu * (1:N)) + sin(3*omega_nu * (1:N)))';
U = algo_fft_1D(u);
figure(1); clf, hold on, plot(u); pause,
figure(2); clf, hold on, plot(0:2/N:1, abs(U(1:end/2+1)), 'r-'); axis([0 1 0 250]); title('segnale stazionario'); pause,
% segnale NON stazionario: 
u = [sin(omega_nu * (1:N/4)) sin(1.5*omega_nu * (N/4+1:N/2)) sin(2*omega_nu * (N/2+1:3*N/4)) sin(3*omega_nu * (3*N/4+1:N))]';
U = algo_fft_1D(u);
figure(3); clf, hold on, plot(u); pause,
figure(4); clf, hold on, plot(0:2/N:1, abs(U(1:end/2+1)), 'r-'); axis([0 1 0 250]); title('segnale NON stazionario'); pause,
%
%
N = 1024;
omega_nu = (2*pi/N) * 0.7*(N/2);
u = (1 + cos(omega_nu * (1:N))) ./ (1 + 4 .* (1:N));
U = algo_fft_1D(u);
figure(5); clf, hold on, plot(u); title('secondo esempio di segnale NON stazionario'); pause,
figure(6); clf, hold on, plot(0:2/N:1, abs(U(1:end/2+1)), 'r-'); title('secondo esempio di segnale NON stazionario'); pause,




%
% Trasformata di Haar 
%
scalnum = [1/sqrt(2) 1/sqrt(2)];
wavenum = [-1/sqrt(2) 1/sqrt(2)];
N = 1024;
omega_nu = (2*pi/N) * 0.01*(N/2);
u = 20  .*  ((0:N-1)./N).^2  .*  (1 - ((0:N-1)./N)).^4  .*  cos(omega_nu * (0:N-1));
figure(7); clf, hold on, plot(u,'b-'); title('segnale originario'); pause,
a = filter(scalnum,[1],u);
b = filter(wavenum,[1],u);
figure(8); clf, hold on, plot([a(2:2:N) b(2:2:N)],'b-'); title('trasformata di Haar - primo livello'); pause,
a2 = filter(scalnum,[1],a(2:2:N));
b2 = filter(wavenum,[1],a(2:2:N));
figure(9); clf, hold on, plot([a2(2:2:N/2) b2(2:2:N/2) b(2:2:N)],'b-'); title('trasformata di Haar - secondo livello'); pause,
%
% facciamola adesso con una matrice sparsa (NB: la sparsità è dovuta al supporto locale delle wavelets):
%
T1 = sparse(1024,1024);
for i=1:N/2  T1(i,:)=[zeros(1,(i-1)*2) 1/sqrt(2) 1/sqrt(2) zeros(1,(N/2-i)*2)]; end
for i=N/2+1:N  T1(i,:)=[zeros(1,(i-N/2-1)*2) 1/sqrt(2) -1/sqrt(2) zeros(1,(N-i)*2)]; end
figure(11); spy(T1); pause
u_haar = T1 * u';
figure(8); plot(u_haar,'r--'); pause,
u_ric = T1' * u_haar;
figure(7); plot(u_ric,'r--'); pause,
%
% trasformata inversa con i filtri:
%
u_ric_filt = filter(scalnum(end:-1:1),[1], upsample(a(2:2:N),2) );
u_ric_filt = u_ric_filt + filter(wavenum(end:-1:1),[1], upsample(b(2:2:N),2));
figure(7); plot(u_ric_filt,'g:'); pause,



%
% contenuto in frequenza delle wavelets
%
haarscalfun = [1/sqrt(2) 1/sqrt(2) zeros(1,256-2)];
haarwavelet = [1/sqrt(2) -1/sqrt(2) zeros(1,256-2)];
figure(10); clf, hold on,
  plot(abs(fft(haarscalfun)).^2,'b-'); 
  plot(abs(fft(haarwavelet)).^2,'r-');  
  plot(abs(fft(haarwavelet)).^2+abs(fft(haarscalfun)).^2,'g-');  
  legend('spettro delle scaling functions','spettro delle wavelets','conserv. energia')
  