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

Biosignals Workshop in Erfurt

Also in the week of DAGA but between the 21st to the 23rd March the Biosignals Workshop took place in the protestant monastery of St. Augustine in Erfurt. The topic was innovative processing of bioelectric and biomagnetic signals and was organized by the technical VDE committees „Biosignals“ and „Magnetic Methods in Medicine“. The DSS group (Christin and Eric) participated with two ...


Read more ...