convolver function for DST (Matlab)

convolver.m

% Michael Bartnett
% mrb402
% Monday November 21, 2011
% convolver.m
%
% Function convolver
%
% Convolves two audio files, a signal and an impulse response
%
% Usage: convolver(IRfilename, SIGfilename, OUTfile, convMethod[, chunkLength])
%
% IRfilename - The name of the mono wave file designated as the impulse response
%
% SIGfilename - The name of the mono wave file designated as the signal
%
% OUTfile - The mono wave file to which to write the result of the convolution
%
% convMethod - 1: Direct convolution, or 2: Fast convolution
% 
% chunkLength - Optional. The length of each chunk in samples for using teh overlap-add method.
%   Defaults to the length of SIGfile
%
function convolver(IRfilename, SIGfilename, OUTfile, convMethod, chunkLength)
    % Verify number of arguments
    if nargin < 4 || nargin > 5
        error('Bad argcount, arguments: IRfilename, SIGfilename, OUTfile, convMethod[, chunkLength]');
    end
 
    % Read the signal and IR from the filesystem
    sigSamples = wavread(IRfilename);
    irSamples = wavread(SIGfilename);
 
    % Ensure chunkLength has at least a default value
    if nargin == 4
        chunkLength = length(sigSamples);
    end
 
    % Find and calculate useful lengths for sample bufers
    irLength = length(irSamples);
    sigLength = length(sigSamples);
    convLength = sigLength + irLength - 1;
    chunkConvLength = chunkLength + irLength - 1;
    bufferOffset = 0;
 
    % Preallocate memory for the output
    convBuffer = zeros(convLength, 1);
 
    % Branch for convolution method
    switch convMethod
    case 1
        % Verify we don't try to read past the length of the signal
        while bufferOffset < sigLength
            % Multiply each sample in IR against each sample in signal
            for i = 1:irLength
                newSampleSlice = sigSamples(1 + bufferOffset:bufferOffset + chunkLength) .* irSamples(i);
                irOffset = bufferOffset + i - 1;
                convBuffer(1 + irOffset:irOffset + chunkLength) = convBuffer(1 + irOffset:irOffset + chunkLength) + newSampleSlice;
            end
            bufferOffset = bufferOffset + chunkLength;
            if (sigLength - bufferOffset) < chunkLength
                chunkLength = sigLength - bufferOffset;
            end
        end
 
    case 2
        % Preallocate chunk buffer
        chunk = zeros(chunkConvLength, 1);
        
        % Resize IR samples to match length of convolved chunk buffer
        resizedIRSamples = irSamples;
        resizedIRSamples(chunkConvLength) = 0;
        
        % Precalculate fft of IR
        irfft = fft(resizedIRSamples);
        
        % Verify we don't try to read past the length of the signal        
        while bufferOffset < sigLength
            % Pull the next chunk out of the source signal and apply FFT
            chunk(1:chunkLength) = sigSamples(1 + bufferOffset:bufferOffset + chunkLength);
            fftChunk = fft(chunk);
            
            % Take inverse fft of freq-domain IR and signal chunk
            %   element-wise multiplied
            backInTime = ifft(irfft .* fftChunk);
            convBuffer(1 + bufferOffset:bufferOffset + chunkConvLength) = convBuffer(1 + bufferOffset:bufferOffset + chunkConvLength) + backInTime;
            
            % Update buffer offset for next chunk
            bufferOffset = bufferOffset + chunkLength;
 
            % If our chunkLength is not an even multiple of the length of the signal, we need to compensate
            if (sigLength - bufferOffset) < chunkLength
                chunkLength = sigLength - bufferOffset;
 
                % If chunk length is zero, we can avoid recalculating all of these
                if chunkLength > 0
                    % Otherwise, the length of the chunk has changed,
                    %   so we need to resize the transformed IR as well
                    chunkConvLength = chunkLength + irLength - 1;
                    chunk = zeros(chunkConvLength, 1);
                    resizedIRSamples = irSamples;
                    resizedIRSamples(chunkConvLength) = 0;
                    irfft = fft(resizedIRSamples);
                end
            end
        end
 
    end
 
    % Normalize the output and write it out
    convBuffer = convBuffer / max(abs(convBuffer));
    wavwrite(convBuffer, 44100, OUTfile);
end