act_wavelet='bior2.2'
decomp_levels=9;
nulovat=90
load Bumps;
insig=Bumps;

inlen=length(insig)
zerobegin=ceil((1-nulovat/100)*inlen);

%spravime DWT
[spec,l] = wavedec(insig,decomp_levels,act_wavelet);

%zistime tresholdovu hranicu odpovedajucu danemu poctu koeficientov
sortspec=sort(abs(spec));
tresh_idx=1+ceil((nulovat/100)*(inlen-1))
treshval=sortspec(tresh_idx)

%vynulujeme zelany pocet koeficientov pomocou hardtresholdingu
zeroed=wthresh(spec,'h',treshval);

%inverzna transformacia
aprox=waverec(zeroed,l,act_wavelet);

subplot(4,1,1); plot(insig);  Ylabel('Vstupny signal');

subplot(4,1,2); plot(spec);   Ylabel('Spektrum');

mytitle= ['Nulované spektrum(',int2str(nulovat),'%)'];
subplot(4,1,3); plot(zeroed); Ylabel(mytitle)

subplot(4,1,4); plot(aprox);  Ylabel('Rekonštruovaný signál');

%vypocitame MSE
MSE=sum((insig-aprox)^2)/inlen