The following matlab function compresses a given vector - zeroing out the terms falling below a user specified (percentage) threshold. This routine is used by all fft and wavelet compression schemes that follow.
function wc=compress(w,r)
% Input is the array w and r, which is a number
% between 0 and 1.
% Output is the array wc where smallest 100r% of the
% terms in w are set to zero (e.g r=.75 means the smallest 75% of
% the terms are set to zero
if (r<0) | (r>1)
error('r should be between 0 and 1')
end;
N=length(w); Nr=floor(N*r);
ww=sort(abs(w));
tol=abs(ww(Nr));
for j=1:N
if (abs(w(j)) < tol)
w(j)=0;
end;
end;
wc=w;
Filtering with FFT. The following Matlab commands filter a given signal. The key Matlab commands are fft and ifft, which compute the fast Fourier transform and its inverse, respectively. Here, the given signal is generated by exponentials, sines and cosines.
>> t=linspace(0,2*pi,2^8); % discretizes [0, 2pi] into 256 nodes >> y=exp(-(cos(t).^2)).*(sin(2*t)+2*cos(4*t)+0.4*sin(t).*sin(50*t)); >> plot(t,y) % generates the graph of the original signal >> fy=fft(y); % computes fft of y >> filterfy=[fy(1:6) zeros(1,2^8-12) fy(2^8-5:2^8)]; % sets fft coefficients >> % to zero for $7 \leq k \leq 256$ >> filtery=ifft(filterfy); % computes inverse fft of the filtered fft >> plot(t, filtery) % generates the plot of the compressed signal
Compression with the FFT. The following routine, called fftcomp.m, uses the FFT and the previous compress.m routine to compress a given signal. The compression rate is part of the input. The routine outputs the graph of the signal, the graph of the compressed signal and the relative l^2 error.
function error=fftcomp(t,y,r)
% Input is an array y, which represents a digitized signal
% associated with the discrete time array t.
% Also input r which is a decimal
% number between 0 and 1 representing the compression rate
% e.g. 80 percent would be r=0.8.
% Outputs are the graph of y and its compression as well as
% the relative error. This routine uses compress.m
%
if (r<0) | (r>1)
error('r should be between 0 and 1')
end;
fy=fft(y);
fyc=compress(fy,r);
yc=ifft(fyc);
plot(t,y,t,yc)
error=norm(y-yc,2)/norm(y)
Use of fftcomp.m for and example on compression with
fft.
>> t=linspace(0,2*pi,2^8); >> y=exp(-t.^2/10).*(sin(2*t)+2*cos(4*t)+0.4*sin(t).*sin(10*t)); >> fftcomp(t,y,0.8) % uses fftcomp with compression rate of 80 percent
function error=daubcomp(t,y,n,r)
% Input is an array y, which represents a digitized signal
% associated with the vector t; n=the number of levels
% (so the number of nodes is 2^n = length of t and the length of y).
% Also input r which is a decimal
% number between 0 and 1 representing the compression rate
% e.g. 80 percent would be r=0.8.
% Output is the graphs of y and its compression, as well as
% the relative error. This routine uses compress.m
% and the Daubechies - 4 wavelets.
%
if (r<0) | (r>1)
error('r should be between 0 and 1')
end;
[c,l]=wavedec(y,n,'db2'); % Matlab's wave decomposition routine
cc=compress(c,r); % compress the signal (compress.m given above)
yc=waverec(cc,l,'db2'); % Matlab's wave reconstruction routine
plot(t,y,t,yc) % plot of the signal and compressed signal
error=norm(y-yc,2)/norm(y) % relative l^2 error
Matlab commands for an Example.
>> t=linspace(0,1,2^8); % discretizes [0,1] into 256 nodes >> y=sin(2*pi*t)+cos(4*pi*t)+sin(8*pi*t)+4*64*(t-1/3).*exp(-((t-1/3)*64).^2)+ 512*(t-2/3).*exp(-((t-2/3)*128).^2); >> daubcomp(t,y,8,0.8)