%% Example 3.15
% This uses the same signal as Examples 3.8 and 3.14. I sample the analog signal f(t) at 64 points where
% f\left(t\right)=16t^2 \text{  }for\text{ }0\text{ }\le t\le \frac{1}{4}, 1 for \frac{1}{4\text{ }}\le t\text{ }\le \frac{1}{2},\text{ }\frac{1}{2}for\text{ }\frac{1}{2}\le t\le \frac{3}{4},and\text{ }2-2t\text{ }for\text{ }\frac{3}{4}\le t\le 1

x = [16*((1:16)/64).^2 ones(1,16) ones(1,16)/2 (2*ones(1,16)-2*(49:64)/64)];
plot(x)        
axis([0 64 -0.2 1.2])
xticks([0:16:64])
title('Signal of length 64 (quadratic, step, linear)')

%
% I compute the discrete Fourier transform of x, replace the 48 smallest components by 0 and take the inverse discrete Fourier transform of the result to obtain a compressed version xcomp. I use sort to sort the components of |y|\text{ }from smallest to largest. The second part of the output gives the corresponding index list, so I can set the 48 smallest elements to 0. I then compute the inverse transform and plot the real part of the result.

y = fft(x);
[A,I] = sort(abs(y));
for j=1:48
    i = I(j);
    y(i)=0;
end
xcomp = ifft(y);
plot(real(xcomp))
axis([0 64 0 1.2])
xticks(0:16:64)
yticks(0:0.5:1.5)
title('Fourier compression')

%
% The relative energy error is about 0.35%.

norm(x-real(xcomp))^2/norm(x)^2

%
% Next I do the same with the  three-scale Daub4 transform. I call the transform w and the compressed transform wcomp.

T6 = daub4mat(64); T5 = [daub4mat(32) zeros(32); zeros(32) eye(32)]; 
T4 = [daub4mat(16) zeros(16,48); zeros(16) eye(16) zeros(16,32); zeros(32) eye(32)];
W3 = T4*T5*T6;
W3s = W3';
w = W3*x';
[B,J] = sort(abs(w));
wcomp = w;
for j = 1:48
    i = J(j);
    wcomp(i)=0;
end
wcomp;
Xcomp=W3s*wcomp;
plot(Xcomp)
axis([0 64 0 1.2])
xticks(0:16:64)
title('Daub4 compression')

%
%
% The relative energy error is about 0.034%, better than the error for Fourier compression by a factor of 10.

norm(x'-Xcomp)^2/norm(x)^2

%

