function compress(x,wvlt,levels,pct)
%
% Example of compression with the DCT, direct calling
% from the the Matlab wavelets toolbox, and Matlab's wdencmp.
%
% To get started, you could run this as follows:
%
% diary filename
% load x.dat
% wvlt='db4';  ????
% levels=13;   ????
% percent=0.2; ????
% compress(x,wvlt,levels,percent)
%
% The program will take perhaps a minute to execute, and generate
% several graphs.
%
% Note that here we use the default dwtmode 'sym'
%
lx=length(x);
normx=norm(x);
fprintf('test signal has length %7.0f\n\n',lx)
figure;
plot(x)
title('signal to be compressed')
fprintf('signal now to be compressed\n\n')
%
% compress with dct
%
fprintf('signal now to be compressed via the dct\n')
fprintf('DCT, amplitudes < %6.3f percent of largest discarded\n',pct)
Xdct=dct(x);
Xdctmax=max(abs(Xdct));
Xthresh=pct*Xdctmax/100;
picked=find(abs(Xdct)<Xthresh);
pdisdct=100*length(picked)/lx;
fprintf('%6.3f percent discarded\n',pdisdct)
Xdctcomp=Xdct;
Xdctcomp(picked)=zeros(size(picked));
xcdct=idct(Xdctcomp);
relnormdct=norm(x-xcdct)/normx;
fprintf('relative error norm %10.3e\n\n',relnormdct)
figure
plot(xcdct)
title('signal compressed with DCT')
%
% compress with Matlab wavelets programs
%
fprintf('compress with Matlab wavelets programs\n')
fprintf('%s wavelet, %3.0f levels, amplitudes < %6.3f percent of largest discarded\n',wvlt,levels,pct)
[W,L]=wavedec(x,levels,wvlt);
Wmax=max(abs(W));
wthresh=pct*Wmax/100;
picked=find(abs(W)<wthresh);
lW=length(W);
pbdisc=100*length(picked)/lW;
fprintf('%6.3f percent discarded from wavelet transform of length %7.0f\n',pbdisc,lW)
Wc=W;
Wc(picked)=zeros(size(picked));
xb=waverec(Wc,L,wvlt);
relnormb=norm(x-xb)/normx;
fprintf('relative error norm %10.3e\n\n',relnormb)
figure
plot(xb)
title('signal compressed with Matlab programs')
%
% compress with Matlab's wdencmp, using above threshold
%
fprintf('compress with Matlabs wdencmp, using above threshold\n')
fprintf('%s wavelet, %3.0f levels\n',wvlt,levels)
keepap=0; % threshold approximation coefficients as well as detail
sorh='h'; % hard thresholding
levdep='gbl'; % threshold the same at all levels
[xc,Cxc,Lxc,perf0,perfl2]=wdencmp(levdep,x,wvlt,levels,wthresh,sorh,keepap);
picked=find(Cxc==0);
pcdisc=100*length(picked)/length(Cxc);
fprintf('%6.3f percent of wavelet transform set = 0\n',pcdisc)
relnormc=norm(x-xc)/normx;
fprintf('relative error norm %10.3e\n',relnormc)
perf0
perfl2
figure
plot(xc)
title('signal compressed with wdencmp')
