% Course in Wavelets, TU Chemnitz, 2012

n=4096;
s=0.01;
r=4;
f=load_signal('Piece-Regular', n);
sigma=s*(max(f)-min(f));
fn=f+sigma*randn(n,1);
% Course in Wavelets, TU Chemnitz, 2012

figure(1)
subplot(2,1,1);
plot(f); axis([1,n,min(f),max(f)]); title('Original');

subplot(2,1,2);
plot(fn); axis([1,n,min(f),max(f)]); title('Noisy');

% Single level DWT
dwtmode('per');
[c,d]=dwt(f,'db4');
[cn,dn]=dwt(fn,'db4');

figure(2)
subplot(2,2,1);
plot(c); axis tight; title('Original, approx. coeffs.');

subplot(2,2,2);
plot(d); axis tight; title('Original, detail coeffs.');

subplot(2,2,3);
plot(cn); axis tight; title('Noisy, approx. coeffs.');

subplot(2,2,4);
plot(dn); axis tight; title('Noisy, detail coeffs.');

% Combined plots, Noisy
figure(3);
[C,L]=wavedec(fn,1,'db4');
subplot(2,1,1);
plot(C); axis tight; title('Noisy, combined coeffs. level 1');

[C,L]=wavedec(fn,3,'db4');
subplot(2,1,2);
plot(C); axis tight; title('Noisy, combined coeffs. level 3');

% Filtering
[C,L]=wavedec(fn,10,'db4');
C(n/2^r+1:end)=0;
ff=waverec(C,L,'db4');

figure(4)
subplot(2,1,1);
plot(f); axis([1,n,min(f),max(f)]); title('Original');

subplot(2,1,2);
plot(ff); axis([1,n,min(f),max(f)]); title('Filtered');

% Thresholding
[C,L]=wavedec(fn,10,'db4');
C(abs(C)<s*max(abs(C)))=0;
ft=waverec(C,L,'db4');

figure(5)
subplot(2,1,1);
plot(f); axis([1,n,min(f),max(f)]); title('Original');

subplot(2,1,2);
plot(ft); axis([1,n,min(f),max(f)]); title('Thresholded');
