%
% compressione (lossy) di segnali audio
%
N=2^16;
u = [hamming(N).*sin(2*pi*(1/164)*(0:N-1)) hamming(N).*sin(2*pi*(1/140)*(0:N-1))];
N = length(u)
en_u = sum(u.^2);  % energia della sequenza "u"
figure(3); plot(u),
wavplay([u' u'], 44100); 
pause
% compressione con wavelets di Haar:
compress = 1;  % 0=no compressione  1=compressione
totflutt = [];
trend = u;
scalnum = [1/sqrt(2) 1/sqrt(2)];
wavenum = [1/sqrt(2) -1/sqrt(2)];
for i=1:5
  utemp = trend;
  trend = filter(scalnum,[1],utemp);
  trend = trend(2:2:end);
  flutt = filter(wavenum,[1],utemp);
  flutt = flutt(2:2:end);
  if compress==0
    totflutt = [flutt totflutt];
  else
    % compressione del segnale: azzero tutti i coefficienti delle fluttuazioni
    totflutt = [zeros(1,length(flutt)) totflutt];
  end
  figure(4); clf, hold on, plot([trend totflutt]); title(['Haar, livello ' num2str(i)]);
  Fc = 44100;
  ttemp = trend;
  % eseguo la trasformata inversa di Haar
  for j=1:i
    ttemp = filter(scalnum(end:-1:1),[1], upsample(ttemp,2) );
    %[sum(N./(2.^(i:-1:i-j+2)))+1 sum(N./(2.^(i:-1:i-j+1)))]
    ttemp = ttemp + filter(wavenum(end:-1:1),[1], upsample(totflutt(sum(N./(2.^(i:-1:i-j+2)))+1:sum(N./(2.^(i:-1:i-j+1)))),2) );
  end
  %
  % calcoliamo in percentuale quanta energia del segnale originario viene mantenuta nelle sequenze "trend"
  % via-via prodotta nella compressione con trasformata di Haar:
  disp(['fraz. energia livello ' num2str(i) ' trasf. Haar: ' num2str(sum(ttemp.^2) / en_u)]);
  %
  wavplay([ttemp' ttemp'], Fc);
  pause
end;
%
% compressione con coiflets  Coif30 :
scalnum = [-.000149963800  .000253561200  .001540245700  -.002941110800  -.007163781900 ...
            .016552066400  .019917804300 -.064997262800  -.036800073600   .298092323500 ...
            .547505429400  .309706849000 -.043866050800  -.074652238900   .029195879500 ...
            .023110777000 -.013973687900 -.006480090000   .004783001400   .001720654700 ...
           -.001175822200 -.000451227000  .000213729800   .000099377600  -.000029232100 ...
           -.000015072000  .000002640800  .000001459300  -.000000118400  -.000000067300]*sqrt(2);
wavenum = [-.000000067300  .000000118400  .000001459300  -.000002640800  -.000015072000 ...
            .000029232100  .000099377600 -.000213729800  -.000451227000   .001175822200 ...
            .001720654700 -.004783001400 -.006480090000   .013973687900   .023110777000 ...
           -.029195879500 -.074652238900  .043866050800   .309706849000  -.547505429400 ...
            .298092323500  .036800073600 -.064997262800  -.019917804300   .016552066400 ...
            .007163781900 -.002941110800 -.001540245700   .000253561200   .000149963800]*sqrt(2);
totflutt = [];
trend = u;
for i=1:5
  utemp = trend;
  trend = filter(scalnum,[1],utemp);
  trend = trend(2:2:end);
  flutt = filter(wavenum,[1],utemp);
  flutt = flutt(2:2:end);
  if compress==0
    totflutt = [flutt totflutt];
  else
    totflutt = [zeros(1,length(flutt)) totflutt];
  end
  figure(4); clf, hold on, plot([trend totflutt]); title(['Coif30, livello ' num2str(i)]);
  Fc = 44100;
  ttemp = trend;
  for j=1:i
    ttemp = filter(scalnum(end:-1:1),[1], upsample(ttemp,2) );
    %[sum(N./(2.^(i:-1:i-j+2)))+1 sum(N./(2.^(i:-1:i-j+1)))]
    ttemp = ttemp + filter(wavenum(end:-1:1),[1], upsample(totflutt(sum(N./(2.^(i:-1:i-j+2)))+1:sum(N./(2.^(i:-1:i-j+1)))),2) );
  end
  %
  % calcoliamo in percentuale quanta energia del segnale originario viene mantenuta nelle sequenze "trend"
  % via-via prodotta nella compressione con trasformata Coif30:
  disp(['fraz. energia livello ' num2str(i) ' trasf. Coif30: ' num2str(sum(ttemp.^2) / en_u)]);
  wavplay([ttemp' ttemp'], Fc);
  pause
end;


%
% rimozione del rumore (de-noising)
%
sigma = 0.2;
n = sigma * randn(1,N);
un = u + n;
soglia = 5*sigma;
figure(3); plot(un), pause,
wavplay([un' un'], 44100); 
totflutt = [];
trend = un;
scalnum = [1/sqrt(2) 1/sqrt(2)];
wavenum = [1/sqrt(2) -1/sqrt(2)];
for i=1:5
  utemp = trend;
  trend = filter(scalnum,[1],utemp);
  trend = trend(2:2:end);
%  trend = trend .* (abs(trend) > soglia);
  flutt = filter(wavenum,[1],utemp);
  flutt = flutt(2:2:end);
  flutt = flutt .* (abs(flutt) > soglia);
  totflutt = [flutt totflutt];
  figure(5); clf, hold on, plot([trend totflutt]); title(['Haar, livello ' num2str(i)]);  pause,
  Fc = 44100;
  ttemp = trend;
  % eseguo la trasformata inversa di Haar
  for j=1:i
    ttemp = filter(scalnum(end:-1:1),[1], upsample(ttemp,2) );
    %[sum(N./(2.^(i:-1:i-j+2)))+1 sum(N./(2.^(i:-1:i-j+1)))]
    ttemp = ttemp + filter(wavenum(end:-1:1),[1], upsample(totflutt(sum(N./(2.^(i:-1:i-j+2)))+1:sum(N./(2.^(i:-1:i-j+1)))),2) );
  end
  wavplay([ttemp' ttemp'], Fc);
  eval(['ttempHaar' num2str(i) ' = ttemp;']);
%  figure(6); clf, hold on,
%             plot(un(4000+[0:400]),'b-');
%             plot(ttemp(4000+[0:400]),'r-');
%  pause,
end;
%
% compressione con coiflets  Coif30 :
scalnum = [-.000149963800  .000253561200  .001540245700  -.002941110800  -.007163781900 ...
            .016552066400  .019917804300 -.064997262800  -.036800073600   .298092323500 ...
            .547505429400  .309706849000 -.043866050800  -.074652238900   .029195879500 ...
            .023110777000 -.013973687900 -.006480090000   .004783001400   .001720654700 ...
           -.001175822200 -.000451227000  .000213729800   .000099377600  -.000029232100 ...
           -.000015072000  .000002640800  .000001459300  -.000000118400  -.000000067300]*sqrt(2);
wavenum = [-.000000067300  .000000118400  .000001459300  -.000002640800  -.000015072000 ...
            .000029232100  .000099377600 -.000213729800  -.000451227000   .001175822200 ...
            .001720654700 -.004783001400 -.006480090000   .013973687900   .023110777000 ...
           -.029195879500 -.074652238900  .043866050800   .309706849000  -.547505429400 ...
            .298092323500  .036800073600 -.064997262800  -.019917804300   .016552066400 ...
            .007163781900 -.002941110800 -.001540245700   .000253561200   .000149963800]*sqrt(2);
totflutt = [];
trend = un;
for i=1:5
  utemp = trend;
  trend = filter(scalnum,[1],utemp);
  trend = trend(2:2:end);
%  trend = trend .* (abs(trend) > soglia);
  flutt = filter(wavenum,[1],utemp);
  flutt = flutt(2:2:end);
  flutt = flutt .* (abs(flutt) > soglia);
  totflutt = [flutt totflutt];
  figure(4); clf, hold on, plot([trend totflutt]); title(['Coif30, livello ' num2str(i)]);  pause,
  Fc = 44100;
  ttemp = trend;
  for j=1:i
    ttemp = filter(scalnum(end:-1:1),[1], upsample(ttemp,2) );
    %[sum(N./(2.^(i:-1:i-j+2)))+1 sum(N./(2.^(i:-1:i-j+1)))]
    ttemp = ttemp + filter(wavenum(end:-1:1),[1], upsample(totflutt(sum(N./(2.^(i:-1:i-j+2)))+1:sum(N./(2.^(i:-1:i-j+1)))),2) );
  end
  wavplay([ttemp' ttemp'], Fc);
  figure(6); clf, hold on,
             lint=150;
             plot(un(4000+[0:lint]),'b-');
             plot(u(4000+[0:lint]),'c-');
             plot(ttemp(4000+[0:lint]),'r-');
             eval(['plot(ttempHaar' num2str(i) '(4000+[0:lint]),''g-'');']);
  pause,
  figure(7); clf, hold on,
             lint=150;
             cf_un = fft(un(4000+[0:lint])); plot(20*log10(abs(cf_un(1:end/2))),'b-');
             cf_u = fft(u(4000+[0:lint])); plot(20*log10(abs(cf_u(1:end/2))),'c-');
             cf_ttemp = fft(ttemp(4000+[0:lint])); plot(20*log10(abs(cf_ttemp(1:end/2))),'r-');
             eval(['cf_ttempHaar = fft(ttempHaar' num2str(i) '(4000+[0:lint])); plot(20*log10(abs(cf_ttempHaar(1:end/2))),''g-'');']);
  pause,
end;
