Discrete Wavelet Transform - Chain processing (Linear Convolution)
For more details, please go to:
http://www.dsprelated.com/showarticle/126.php
Important, to run the DWT chain processing version you also need:
recorreder.m - http://www.dsprelated.com/showcode/43.php
formfilters.m (chain version) - http://www.dsprelated.com/showcode/44.php
formfilterdwt.m (chain version) - http://www.dsprelated.com/showcode/45.php
getsincdwt.m (chain version) - http://www.dsprelated.com/showcode/46.php
% ----------------------------------------------------
% Title: Discrete Wavelet Transform
% Author: David Valencia
% UPIITA IPN 2010
%
% Posted in DSPrelated.com
% http://www.dsprelated.com/showcode/47.php
%
% Computes the Discrete Wavelet Transform of
% n levels (and its branches) from a n sample input.
% Base filters can be of any lenght.
%
% Generalized DWT filter bank computed via the
% chain processing method (linear convolution)
%
% Dependencies:
% formfilters.m - http://www.dsprelated.com/showcode/44.php
% formfiltersdwt.m - http://www.dsprelated.com/showcode/45.php
% upsample2.m - http://www.dsprelated.com/showcode/10.php
% recorreder.m - http://www.dsprelated.com/showcode/43.php
% getsincdwt.m - http://www.dsprelated.com/showcode/46.php
%
% Revisions:
% v1.0a Commented and translated in English
% ----------------------------------------------------
close all; clear all; clc;
% ######################################
disp('== CHAIN PROCESSING DWT ==')
% Filter parameters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db2';
if(typeofbasis == 'b')
[Rf,Df] = biorwavf(typbior);
[h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
[h0,h1,g0,g1] = wfilters(typor);
end;
% Input values, change as needed
N = 16;
x = (1:N);
L = length(h0);
figure(1);
subplot(2,1,1);
stem(x);
xlabel('Original signal');
% Filter bank parameters, change as needed
n_stages = 3; %Of the low-pass part
% Checked for the minimal case: DWT = DWPT n_stages = 1
n_branches = n_stages + 1; %Of the low-pass part
niter = 300; %Number of iterations
% Define input buffer dimensions
hx = formfilters(n_stages, 1, h0, h1);
lhx = length(hx);
% Input buffer vector
inputbuffer = zeros(lhx,1);
% Synthesis input buffer
sbuffer = zeros(n_branches,lhx);
% This vector will store the length of each branch filter
Lfiltros = zeros(n_branches, 1);
% This vector will store the decimate factor for each branch
dec_factors = zeros(n_branches, 1);
% This vector will store the synchronization factor for each branch
sinc_factors = zeros(n_branches, 1);
% Cycle to get the synchronization factors
for j = 0:n_branches-1
try
sinc_factors(j+1) = getsincdwt(n_stages, j+1, L);
catch error
disp('ERROR: Error, experimentation is needed for this case');
disp('The results may not be correct');
% return;
end
end
%% Filter matrices formation
% This matrix will store each of the branch filters as vectors
H = zeros(n_branches, lhx); %Malloc
G = zeros(n_branches, lhx); %Malloc
for i = 0:n_stages
if i >= n_stages-1
dec_factors(i+1) = 2^n_stages;
else
dec_factors(i+1) = 2^(i+1);
end
hx = fliplr(formfiltersdwt(n_stages,n_stages+1-i,h0,h1));
gx = fliplr(formfiltersdwt(n_stages,n_stages+1-i,g0,g1));
lhx = length(hx);
lgx = length(gx);
Lfiltros(i+1) = lhx;
H(i+1,1:lhx) = hx;
G(i+1,1:lhx) = gx;
end
% Debug code
disp('Dimensiones de los filtros')
Lfiltros
disp('Factores de diezmado')
dec_factors
disp('Matriz de filtros análisis: ');
H
disp('Matriz de filtros sÃntesis: ');
G
chanbufftestigo = zeros(n_branches,1);
outputbuffer = zeros(1,niter);
%% MAIN CYCLE
for i = 1:niter %General FOR (for iterations)
i % Print iteration number
% Shifting of input buffer (sequentially)
inputbuffer = recorreder(inputbuffer,1);
if i<=N
inputbuffer(1)=x(i);
else
inputbuffer(1)=0;
end
inputbuffer %Print buffer (debug)
%% Analyisis stage (h0 and h1)
% The following cycle will select each of the branches for processing
% If the iteration counter is a multiply of the decimation factor,
% the convolution is calculated, otherwise, zeros are sent to the
% channel
y = zeros(n_branches,1); %Clear the output buffer
for j = 0:n_branches-1
fprintf('\nBranch: %d, Decimating factor: %d\n',j+1,dec_factors(j+1));
fprintf('\nFilter length: %d\n',Lfiltros(j+1));
if mod(i,dec_factors(j+1)) == 1
fprintf('i = %d is a multiply of the factor: %d\n',i,dec_factors(j+1));
tempfilt = H(j+1,1:Lfiltros(j+1));
% If the current iteration (clock cycle) is a multiply of the
% factor, the convolution is computed
for k=1:Lfiltros(j+1) %convolution
y(j+1,1) = y(j+1,1) + inputbuffer(k)*tempfilt(k);
end;
end
disp('Results in the channel');
chanbuff(:,1) = y %Debug printing
end
%% Synthesis stage (g0 and g1)
% Shifting and interpolation of the synthesis buffer
for j = 1:n_branches
sbuffer(j,:) = recorreder(sbuffer(j,:),1);
% Interpolation of the synthesis stage inputs
if mod(i,dec_factors(j)) == 1
fprintf('Inserting sample to the synthesis branch %d with dec_factor %d\n',j,dec_factors(j));
sbuffer(j,1) = chanbuff(j,1);
else
fprintf('Inserting a zero to the synthesis branch %d with dec_factor %d\n',j,dec_factors(j));
sbuffer(j,1) = 0;
end
end
disp('Synthesis buffer matrix');
sbuffer
xnbuff = zeros(n_branches,1);
for j=0:n_branches-1
fprintf('Branch: %d , dec_factor = %d',j+1,dec_factors(j+1));
% === BRANCH SYNCHRONIZATION ===
% The branches (except the two last ones down the bank filter)
% will only take a part of the vector that corresponds with their
% channel, taking 'Lx' elements, being Lx the number of coefficient
% that their filter has. An offset is applied
if j < n_branches - 2
[m,n] = size(sbuffer);
tempsinth = sbuffer(j+1,n-Lfiltros(j+1)+1:n) %TESTING
else % The lowest branches take the whole vector
tempsinth = sbuffer(j+1,1+sinc_factors(j+1):Lfiltros(j+1)+sinc_factors(j+1))
end
for k = 1 : Lfiltros(j+1) %Convolution
xnbuff(j+1,1) = xnbuff(j+1,1) + tempsinth(k)*G(j+1,k);
end
end
xnbuff
outputbuffer(i) = sum(xnbuff);
disp('Output buffer: ');
if i > N
disp(outputbuffer(1,i-N+1:i))
figure(1);
subplot(2,1,2);
stem(outputbuffer(1,i-N+1:i));
xlabel('Recovered signal');
else
disp(outputbuffer(1,1:N))
figure(1);
subplot(2,1,2);
stem(outputbuffer(1,1:N));
xlabel('Recovered signal');
end
pause;
end
% END OF PROGRAM >>>
