Discrete Wavelet Transform - Chain processing (Linear Convolution)

David Valencia December 5, 20101 comment Coded in Matlab

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 >>>

Comments:

David VP
Said:
The links to the other neccesary code snippets have been corrected now. Thanks for the report.
8 years ago
0
Reply
Sorry, you need javascript enabled to post any comments.
Sorry, you need javascript enabled to post any comments.