Fast Convolution without Additional Delay

 

written by Gerhard Schmidt, Anton Namenas, and Seedo Eldho Paul

 

The section is currently in progress ...

Please enter the password to see the content

 

Code Example

As usual you can find below some matlab code that shows a comparison between a pure time-domain convolution and the described mixed-domain convolution.

Convolutions examples

%**************************************************************************
% Basic parameters
%**************************************************************************
N     = 588;                % Filter lenght
r     =  64;                % Frameshift
N_FFT = 128;                % FFT size

%**************************************************************************
% Input signal (white Gaussian noise)
%**************************************************************************
x = randn(5000,1);          % Input signal

%**************************************************************************
% Impulse response (white Gaussian noise)
%**************************************************************************
h = randn(N,1);             % Impulse response

%**************************************************************************
% Pure time-domain convolution
%**************************************************************************
y = conv(x,h);              % output signal

%**************************************************************************
% Mixed-domain convolution
%**************************************************************************

    %**********************************************************************
    % Initialization
    %**********************************************************************
    x_td_buffer     = zeros(N_FFT,1);
    h_td            = h(1:N_FFT);
    y_td            = zeros(size(x));
    y_fd_res_buffer = zeros(size(x));
    y_td_curr       = zeros(N_FFT,1);
    k_fd            = 0;
    M               = ceil((N-2*r)/r);
    X_fd_buffer     = zeros(N_FFT/2+1,M);
    H_fd_buffer     = zeros(N_FFT/2+1,M);
    
    %**********************************************************************
    % Fill the frequency-domain filter coefficients
    %**********************************************************************
    h_ind = N_FFT;
    
    % Loop over all frames
    for m = 1:M
        
        %******************************************************************
        % Reset impulse response vector
        %******************************************************************
        h_curr = zeros(r,1);
        
        %******************************************************************
        % Fill the vector with the corresponding parts of the impulse res.
        %******************************************************************
        for n = 1:r
            h_ind = h_ind + 1;
            if (h_ind <= N)
                h_curr(n) = h(h_ind);
            end;
        end;
        
        %******************************************************************
        % Compute FFT
        %******************************************************************
        H_curr           = fft(h_curr,N_FFT);
        H_fd_buffer(:,m) = H_curr(1:N_FFT/2+1);
    end;
    
    %**********************************************************************
    % Main loop
    %**********************************************************************
    for k = 1:length(x)
        
        %******************************************************************
        % Update counters
        %******************************************************************
        k_fd = k_fd + 1;
        if (k_fd > r)
            k_fd = 1;
        end;
        
        %******************************************************************
        % Update the buffer
        %******************************************************************
        x_td_buffer(1:end-1) = x_td_buffer(2:end);
        x_td_buffer(end)     = x(k);
        
        %******************************************************************
        % Generate time-domain based part of the output 
        %******************************************************************
        y_td(k) = x_td_buffer(end:-1:1)' * h_td + y_fd_res_buffer(k_fd);
        
        %******************************************************************
        % Start subsampled processing
        %******************************************************************
        if (k_fd == r)
            
           %***************************************************************
           % Save the result of the previous background processing
           %***************************************************************
           y_fd_res_buffer = y_td_curr(r+1:end);
            
           %***************************************************************
           % Compute FFT on the input, if one full frame is available
           %*************************************************************** 
           X_curr = fft(x_td_buffer,N_FFT);
           X_curr = X_curr(1:N_FFT/2+1);
           
           %***************************************************************
           % Update the FFT buffer
           %***************************************************************
           X_fd_buffer(:,2:M) = X_fd_buffer(:,1:M-1);
           X_fd_buffer(:,1)   = X_curr;
     
           %***************************************************************
           % Compute the frequency-domain convolution output
           %***************************************************************
           Y_fd = zeros(N_FFT/2+1,1);
           for m = 1:M
               Y_fd = Y_fd + X_fd_buffer(:,m) .* H_fd_buffer(:,m); 
           end;
           
           %***************************************************************
           % Compute inverse FFT of the output
           %***************************************************************
           Y_fd_curr       = [Y_fd; conj(Y_fd(end-1:-1:2))];
           y_td_curr       = ifft(Y_fd_curr);
           
        end;      
    end;
   
%**************************************************************************
% Show the result of both convolutions
%**************************************************************************
offset = 2;
lw     = 1;

figure(1);
plot(y(1:500),'b','LineWidth',lw);
hold on
plot(y_td(1:500)+offset,'r','LineWidth',lw);
grid on
hold off
legend('Output (time domain)', ...
       ['Output (mixed domain) + ',num2str(offset)]);
xlabel('Samples')


 

Basic Idea of RED

RED (robust and efficient DSP) is a collection of algorithms that should help engineers to implement basic signal processign tasks in - as the name states - a robust and efficient way. This includes mainly numerical aspects.

 

 

Here you can find the chapter on fast convolution without additional delay as a single pdf file. If you would like to download all chapters please go to the RED overview page and get the pdf document that contain all chapters from there.

 

Authors of this Section:

Anton Namenas
Seedo Eldho Paul
Gerhard Schmidt

 

Related Sections/Chapters:

 

Related Downloads:

Pdf file of this section
   
Matlab script for a convolution without additional delay using partly FFT-based approaches.

Contact

Prof. Dr.-Ing. Gerhard Schmidt

E-Mail: gus@tf.uni-kiel.de

Christian-Albrechts-Universität zu Kiel
Faculty of Engineering
Institute for Electrical Engineering and Information Engineering
Digital Signal Processing and System Theory

Kaiserstr. 2
24143 Kiel, Germany

Recent News

New PhDs in the DSS Team

Since January this year we have two new PhD students in the team: Elke Warmerdam and Finn Spitz.

Elke is from Amsterdam and she works in the neurology department in the university hospital in the group of Prof. Maetzler. Her research topic is movement analysis of patients with neurologic disorders. Elke cooperates with us in signal processing related aspects of her research. Elke plays ...


Read more ...