% fig31
clear;clc;close all;
N=512;
err1=[];
err2=[];

z=zeros(1,N);
z(1,[129:256])=sin((abs([128:255]-128).^1.7)/128);
z(1,[385:448])=sin((abs([384:447]-128).^2)/128);
subplot(221);plot(z);axis([0 512 -1.1 1.1]);
title('Original Signal');

for K=1:400,
% Fourier compression
z_hat=fft(z);
zz=fliplr(sort(abs(z_hat)));

thres=zz(1,K);

for k=1:N,
   if (abs(z_hat(k))<thres)
      z_hat(k)=0;
   end;   
end;   

w1=ifft(z_hat);


% Wavelet compression (db6)
[c,l] = wavedec(z,4,'db6');
zz=fliplr(sort(abs(c)));

thres=zz(1,K);

for k=1:N,
   if (abs(c(k))<thres)
      c(k)=0;
   end;   
end;

w2 = waverec(c,l,'db6');

% Relative error
err1=[err1 norm(z-w1)./norm(z)];
err2=[err2 norm(z-w2)./norm(z)];

% Plotting
if K==75
subplot(222);plot(real(w1));axis([0 512 -1.1 1.1]);
title('Fourier compression,K=75');
subplot(223);plot(w2);axis([0 512 -1.1 1.1]);
title('D6 wavelet compression,K=75');
end;
end;
x=length(err1);
subplot(224);plot(1:x,err1,':',1:x,err2);
title('Relative error');
legend('Fourier basis','D6 basis');