clc, close all, clear all;

%% generate a signal
fx    = 10e3;
fc    = 25*fx;
Tstop = 20/fx;
Amp   = 1;
t     = 0:(1/fc):Tstop;
rand  = 0.25 * randn(size(t));                                             % additive white noise

signal = Amp * sin(2*pi*fx*t) + rand;

% figure()
%    plot(t, signal);
%    title('segnale - dominio del tempo');
%    xlabel('time [s]');
%    ylabel('amplitude');

%% decompose signal using DWT (discrete wavelet transform)
level  = 5;
wname  = 'db30';
type   = 'a';
[C, L] = wavedec(signal, level, wname);                                    % performs a multilevel 1-D wavelet analysis. C contains the wavelet decomposition. L contains the number of coefficients by level
xN     = waverec(C, L, wname);                                             % performs a multilevel 1-D wavelet reconstruction

cAppr  = appcoef(C, L, 'db30', level);                                          % calcola l'ultima approssimazione
[cD1 cD2 cD3 cD4] = detcoef(C, L, [1,2,3,4]);                                   % calcola i dettagli

A      = wrcoef(type, C, L, wname, level);                                 % reconstructs the coefficients of a 1-D signal
                                                                           % type determines whether approximation ('a') or detail ('d') coefficients are reconstructed.
figure(),
   for i= 1:level
       DD(i,:) = wrcoef('d', C, L, wname, i);
       subplot(level+2, 1, i);
       plot(DD(i,:));
       stringa = sprintf('dettaglio D%d', i);
       title(stringa);
   end

   subplot(level+2, 1, level+1), plot(A);
   stringa = sprintf('approssimazione A%d', level);
   title(stringa);
   subplot(level+2, 1, level+2), plot(signal, 'r');
   title('segnale');

%% denoise the signal
signald = wden(signal, 'rigrsure', 's', 'sln', level, wname);              % automatic 1-D denoising using wavelets
                                                                           % signal - threshold technique - soft thresholding - threshold rescaling using single estimate of noise based on the first level coefficent - decomposition level - wavelet
figure(),
   subplot(2, 1, 1);
   plot(signal);
   axis tight;
   grid on;
   title('noisy signal');

   subplot(2, 1, 2);
   plot(signald);
   axis tight;
   grid on;
   title('denoised signal');

%%
DN = detcoef(C, L, 1);                                                     % extracts the detail coefficients at level N from the wavelet decomposition structure [C,L]