APPENDIX C
MATLAB ROUTINES
C.1 GENERAL COMPRESSION ROUTINE
The following Matlab function compresses a given vector, thereby 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;
C.2 USE OF MATLAB’S FFT ROUTINE FOR FILTERING AND COMPRESSION
Filtering with FFT. The following Matlab commands are needed for Example 3.6, which discusses filtering the high-frequency components from a given signal. The key Matlab commands are fft and ifft, which compute the fast Fourier transform and its inverse, respectively.
>> 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 l2 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 Example 3.7 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
C.3 SAMPLE ROUTINES USING MATLAB’S WAVELET TOOLBOX
Matlab Commands Needed for Compression and Filtering with Wavelets. The following Matlab routine decomposes a signal into a wavelet decomposition using Matlab’s wavelet package. The routine compresses this decomposition and then reconstructs the compressed signal. This routine used two key Matlab commands, wavedec and waverec, for the decomposition and reconstruction. Different wavelets can be used with this routine. We use Daubechies 4-coefficient wavelets (indicated by the parameter db2). Higher-order wavelets can be used (denoted dbn, where n is an integer from 1 to 50; n = 1 denotes Haar wavelets). Inputs for this routine are the signal, y, the associated time nodes, t, the number of levels, n, in the discretization, and the compression rate r. The routine outputs the graphs of the signal and the compressed signal as well as the relative l2 error.
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 Figure 5.14
>> 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)
The same routine with db2 replaced by db1 (for Haar) is used for Example 4.15.
C.4 MATLAB CODE FOR THE ALGORITHMS IN SECTION 5.2
The following Matlab routine (called dec) will take a given signal as input and returns the plot of the Vj-component of the signal where j is prescribed by the user. The wavelet decomposition down to level j is also returned. This is not intended to be a professional code (the Matlab wavelet toolbox provides professional code). Rather, the intent here is to show how Matlab can be used to encode the decomposition and reconstruction algorithms given in Section 5.2.
Decomposition
function w=dec(f,p,NJ, Jstop)
%Inputs: f = data whose length is 2^NJ, where NJ=number of scales.
% p = scaling coefficients
% Jstop = stopping scale; program will decompose down
% to scale level Jstop.
%Outputs: w=wavelet coefficients down to level W-Jstop
% the first 1:2^Jstop entries of w is the V-Jstop projection
% of f. The rest are the wavelet coefficients.
N=length(f); N1=2^NJ;
if ˜(N==N1)
error(’Length of f should be 2^NJ’)
end;
if (Jstop <1)|(Jstop>NJ)
error(’Jstop must be at least 1 and <= NJ’)
end;
L=length(p);
pf=fliplr(p);
q=p; q(2:2:L) = -q(2:2:L);
a=f;
t=[];
for j=NJ:-1:Jstop+1
n=length(a);
a=[a(mod((-L+1:-1),n)+1) a]; % make the data periodic
b=conv(a,q); b=b(L+1:2:L+n-1)/2;
a=conv(a,pf); a=a(L:L+n-1)/2; % convolve
ab=a(1:L); a=[a(L+1:n) ab]; % periodize
a=a(2:2:n); % then down-sample
t=[b,t];
end;
w=[a,t]; JJ=2^(Jstop);
ww=[w(JJ) w(1:JJ)]; % returns a periodic graph
tt=linspace(0,1,JJ+1);
if tt=2 % for Haar, the following plot routine returns a block graph
ll=length(tt);
ta=[tt; tt]; tt=ta(1:2*ll);
wa=[ww; ww]; ww=wa(1:2*ll);
ww=[ww(2*ll) ww(1:2*ll-1)];
end;
plot(tt,ww)
Here is the Matlab session that generates Figure 5.11 on the V4 component of a given signal.
>> t=linspace(0,1,2^8); % discretize the unit interval into 2^8 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); % Sample signal >> p=[0.6830 1.1830 0.3170 -0.1830] % Coefficients for Daubechies -4. >> w=dec(y,p,8,4); % decomposes the signal y from level 8 down to level 4
Reconstruction. The following code takes the wavelet decomposition of a given signal down to level j (where j is prescribed by the user) and reconstructs the signal to top level.
function y=recon(w,p,NJ, Jstart)
%Inputs: w = wavelet coefficients length is 2^NJ, where NJ=number of scales.
% ordered by resolution (from lowest to highest).
% p = scaling coefficients
% Jstart = starting scale; program will reconstruct
% starting with V_Jstart and ending with NJ
%
%Outputs: y=reconstructed signal at V_NJ with a corresponding plot
%
N=length(w); Nj=(2^Jstart);
if ˜(N==2^NJ)
error(’Length of w should be 2^NJ’)
end;
if (Jstart <1)|(Jstart>NJ)
error(’Jstop must be at least 1 and <= NJ’)
end;
L=length(p);
q=fliplr(p);
a=w(1:Nj);
for j=Jstart:(NJ-1)
b=w(Nj+1:2*Nj);
m=mod((0:L/2-1),Nj)+1;
Nj=2*Nj;
ua(2:2:Nj+L)=[a a(1,m)]; % periodize the data and upsample
ub(2:2:Nj+L)=[b b(1,m)]; % periodize the data and upsample
ca=conv(ua,p); ca=[ca(Nj:Nj+L-1) ca(L:Nj-1)]; % convolve with p
cb=conv(ub,q); cb=cb(L:Nj+L-1); % convolve with q
cb(1:2:Nj)=-cb(1:2:Nj); % sign change on the odd entries
a=ca+cb;
end;
y=a;
yy=[y(N) y]; % periodize the data
t=linspace(0,1,N+1);
if L==2 % in the Haar case, return a block-style graph
ll=length(t);
ta=[t; t]; t=ta(1:2*ll);
ya=[yy; yy]; yy=ya(1:2*ll);
yy=[yy(2*ll) yy(1:2*ll-1)];
end;
plot(t,yy)
The following Matlab session compresses the signal f using 80% compression and reproduces Figure 5.14.
>> wc=compress(w,0.8); >> t=linspace(0,1,2^8); % discretize the unit interval into 2^8 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); % Sample signal >> p=[0.6830 1.18300.3170 -0.1830] % Coefficients for Daubechies -4. >> w=dec(y,p,8,1); %decomposes the signal y from level 8 down to level 1 >> wc=compress(w,0.8); %compresses the wavelet coefficients by 80 percent >> %(compress is the routine given at the beginning of this section) >> yc=recon(wc,p,8,1); % reconstructs from level 1 to level 8 >> % from the compressed wavelet coefficients wc