Sign in

Not a member? | Forgot your Password?

Search code

Search tips

Find us on Facebook!





Free PDF Downloads

A Quadrature Signals Tutorial: Complex, But Not Complicated

Understanding the 'Phasing Method' of Single Sideband Demodulation

Complex Digital Signal Processing in Telecommunications

Introduction to Sound Processing

C++ Tutorial

Introduction of C Programming for DSP Applications

Fixed-Point Arithmetic: An Introduction

Cascaded Integrator-Comb (CIC) Filter Introduction

FFT Spectral Analysis Software

See Also

Embedded SystemsFPGA

DSP Code Sharing > Fast convolution

Fast convolution

Language: Matlab

Processor: Not Relevant

Submitted by Hernán Ordiales on Nov 9 2010

Licensed under a Creative Commons Attribution 3.0 Unported License

Fast convolution


 

Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.

Uses FFT to perform the convolution by multiplication in the frequency domain. To reduce in this way, the computation time.

Works with mono or stereo input signals.

 
%
%   Copyright (C) 2004 by Hernán Ordiales
%   <audiocode@uint8.com.ar>
%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%

function fast_conv( file, impulse, output )
% Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
% USE: fast_conv( 'source.wav', 'impulse_response.wav', 'output.wav' );

h_stereo = true;
x_stereo = true;

clip_factor = 1.01; % clip factor default value

% impulse response:
[ h, Fs1 ] = wavread( impulse ); % 2 columns matrix, one per channel
h1 = h(:,1); % col 1 -> left channel
h_channels = size(h); h_channels = h_channels(2);
if( h_channels == 2 )
        h2 = h(:,2); % col 2 -> right channel
else
        h2 = h(:,1); % col 1 -> left channel
        h_stereo = false;
end

% file to process:
[ x, Fs2 ] = wavread( file );
x1 = x(:,1); % col 1 -> left channel
x_channels = size(x); x_channels = x_channels(2);
if( x_channels == 2 )
        x2 = x(:,2); % col 2 -> right channel
else
        x2 = x(:,1); % col 1 -> left channel
        x_stereo = false;
end

% if sample rates are the same
if( Fs1 == Fs2 )

        % searches for the amount of points required to perform the FFT
        L = length(h1) + length(x1) - 1; % linear convolution length

        disp('Processing...');

       
        % Use nextpow if you have it to replace below code -> N = 2^nextpow(L);
        N = 2;
        while( N<L ) % first 2 pow greater than L
                N = N*2;
        end
       

        % Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
        % common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
       

        % FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
        H1 = fft( h1, N ); % Fourier transform of the impulse
        X1 = fft( x1, N ); % Fourier transform of the input signal
       
        F_y1 = H1 .* X1; % spectral multiplication

        y1 = ifft( F_y1 ); % time domain again

        % audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
        y1 = y1/( max(y1)*clip_factor ); % to avoid clipping

        % if one of the files is stereo, process the other channel
        if( h_stereo || x_stereo )
                H2 = fft( h2, N ); % Fourier transform of the impulse
                X2 = fft( x2, N ); % Fourier transform of the input signal

                F_y2 = H2 .* X2; % spectral multiplication
                y2 = ifft( F_y2 ); % time domain again

                y2 = y2/( max(y2)*clip_factor ); % to avoid clipping

                y = [ y1 y2 ]; % stereo output (two channels)
        else
                y = y1; % mono output
        end

        wavwrite( y, Fs1, output );

        disp( 'Convolution success.' );

else

        disp('Error: files has different sample rate.');

end
 
 
Rate this code snippet:
0
Rating: 0 | Votes: 0
 
   
 
posted by Hernán Ordiales



Comments


No comments yet for this code


Add a Comment
You need to login before you can post a comment (best way to prevent spam). ( Not a member? )