Multichannel signal filter with multichannel impulse response (efficiently)

Objective:

Filter 16-channel audio files with corresponding 16-channel impulse response (IR) data (obtained using MATLAB IR measurer app and stored as .mat file in the default format) efficiently. Each channel is filtered independently by one IR only, on a 1-1 (...16-16) basis.

Tried:

Running this in a nested loop works, but is relatively slow for filtering lots of audio files. Simple example:

% load impulse responses (format is according to MATLAB app output)
load("measured_ir_data.mat")

% get audio render filenames (filepath variable is the system path to the
% input audio files)
filelist = string({dir(fullfile(filepath ,"*.wav")).name}.');

% loop over files and filter with each IR
for ii = length(filelist):-1:1
    filename = filelist(ii);
    [inAudio, sampleRate] = audioread(filename);
    
    % loop over channels
    for jj = size(inAudio, 2):-1:1
        filtAudio(:, jj) = filter(measurementData.ImpulseResponse(jj).Amplitude,...
                                  1, inAudio(:, jj), [], 1);
    end
    
    % the filtered audio is then further processed and written to file on
    % each loop iteration
end

Perhaps this could be done more efficiently using a dsp.frequencydomainfirfilter-system-object? But that seems to want to filter each input channel with every filter, which is not what is needed here.

Perhaps there is another, more efficient approach using filtering, that does not involve resorting to matrix convolution?

[Hitesh]

Hi @Michael

Can you share "measured_ir_data.mat" and ".wav" file for reporducing the time delay so that I can reproduce the time delay and explore a more efficient approach?

[Michael]

Apologies, I should have provided a simpler example and included some dummy data to facilitate.

% generate some dummy impulse responses
dummyIR = randn(25*48e3, 16)/10.*exp(-1e-4.*linspace(0, 25*48e3 - 1/48e3, 25*48e3).');

fileNumber = 200; % this is arbitrary, just to indicate the quantity of files involved

% loop over files and filter with each IR
for ii = fileNumber:-1:1
    % generate some dummy audio data on each loop iteration
    dummyData = randn(25*48e3, 16)/10;

    % loop over channels
    for jj = size(dummyData, 2):-1:1
        filtAudio(:, jj) = filter(dummyIR(:, jj),...
                                  1, dummyData(:, jj), [], 1);
    end
    
    % the filtered audio is then further processed and written to file on
    % each loop iteration
end

Answer 1

[jibrahim]

Hi Michael,

Staritng in R2023B, dsp.FrequencyDomainFilter now supports speciying multiple filters. you can leverage that ability to perform what you want. You should see very significant speedups for long impulse responses. For example, we can revisit your code:

% Generate some dummy impulse responses
numChans = 16;
Fs = 48e3;
IRLen = round(5*Fs); % 5 seconds - in samples
dummyIR = randn(IRLen, numChans)/10.*exp(-1e-4.*linspace(0, IRLen - 1/Fs, IRLen).');

% This is arbitrary, just to indicate the quantity of files involved
numFiles = 2; 
audioLen = round(25*Fs); % 25 sec - in samples
filtAudio = zeros(audioLen,numChans,numFiles); % results go here

% generate some dummy audio data on each loop iteration
dummyDataAll = randn(audioLen, numChans,numFiles)/10;

tic;
% loop over files and filter with each IR
for ii = numFiles:-1:1   
    dummyData = dummyDataAll(:,:,ii);
    % loop over channels
    for jj = size(dummyData, 2):-1:1
        filtAudio(:, jj,ii) = filter(dummyIR(:, jj),...
                                  1, dummyData(:, jj), [], 1);
    end
end
t1 = toc % on my machine, 268.7142 seconds

Now do the same, but with dsp.FrequencyDomainFilter.

filtAudio2 = zeros(audioLen,numChans,numFiles);
f = dsp.FrequencyDomainFIRFilter(Numerator = dummyIR.',NumPaths=1,SumFilteredOutputs=false);

tic;
% loop over files and filter with each IR
for ii = numFiles:-1:1
    dummyData = dummyDataAll(:,:,ii);
    filtAudio2(:,:,ii) = squeeze(f(dummyData));
    reset(f)
end
t2 = toc % on my machine, 2.8056 sec

On my machine, I see a 95X speed-up.

When comparing results, account for latency caused by overlap-add of the freuqncy-domain filter. You can control latency by using PartitionForReducedLatency on the object. Partioning will reduce latency at the expense of computational speed. Padding your signals with zeros should allow we to get the entire response while using dsp.FrequencyDomainFIRFilter

% Compare results. Account for latency.
f1 = filtAudio(1:end-f.Latency,:,:);
f2 = filtAudio2(f.Latency+1:end,:,:);
norm(f1(:)-f2(:)) %  1.3176e-11 (good)

If you are working with files on disk, please note that you can simoplify your code and get more speedups by using audioDatastore. Datastores allow you to convientely point to multiple files and process them in parallel if you have access to Parallel Processing Toolbox.

For example, on my machine, I assume I have 20 files uder a subfolder named myfiles:

% Create the datastore
ads = audioDatastore(fullfile(pwd,"myfiles"));

% Create a transform datastore that applies filtering after reading data
% from file:
fds = transform(ads,@(x)myFilter(x,f));

% Process al the files
tic;
Y1 = readall(fds,UseParallel=false);
toc % on my machine, 18.35 s
Y1 = reshape(Y1,[],numChans,numel(ads.Files));

% Now use parallel processin toolbox to read in parallel on different CPU
% workers on your machine
tic;
Y2 = readall(fds,UseParallel=true);
toc % on my machine, 9.8 s

% Results have been concatenated into one big matrix. Reshape and compare.
Y2 = reshape(Y2,[],numChans,numel(ads.Files));
norm(Y1(:)-Y2(:))

function y = myFilter(x,f)
    y = squeeze(f(x));
end

[Michael]

Thanks @jibrahim for your comprehensive answer. I have tried this out, and there seems to be a potential problem: the latency on the filter obtained is very large (for my 4-s IRs, the latency obtained from f.Latency seems to be 4-s). Attempting to use the PartitionForReducedLatency option fails because this option does not (currently) support multiple filters in the filter object.

Does the latency only matter if the intended application is real-time? How would it affect the output signal in a post-processing situation, if at all?

[jibrahim]

(Note that support for multiple filters with reduced latency is supported starting in R2024A)

Hi Michael, you can pad the input with zeros and then remove the leading output section to get the same results as filter:

%Generate some dummy impulse responses
numChans = 16;
Fs = 48e3;
IRLen = round(5*Fs); % 5 seconds - in samples
dummyIR = randn(IRLen, numChans)/10.*exp(-1e-4.*linspace(0, IRLen - 1/Fs, IRLen).');

%This is arbitrary, just to indicate the quantity of files involved
numFiles = 2; 
audioLen = round(10*Fs); % 25 sec - in samples
filtAudio = zeros(audioLen,numChans,numFiles);
% generate some dummy audio data on each loop iteration
dummyDataAll = randn(audioLen, numChans,numFiles)/10;
tic;
% loop over files and filter with each IR
for ii = numFiles:-1:1   
    dummyData = dummyDataAll(:,:,ii);
    % loop over channels
    for jj = size(dummyData, 2):-1:1
        filtAudio(:, jj,ii) = filter(dummyIR(:, jj),...
                                  1, dummyData(:, jj), [], 1);
    end
end
t1 = toc % 84.6570

filtAudio2 = zeros(audioLen,numChans,numFiles);
f = dsp.FrequencyDomainFIRFilter(Numerator = dummyIR.',NumPaths=1,SumFilteredOutputs=false);
tic;

% loop over files and filter with each IR
for ii = numFiles:-1:1
    dummyData = dummyDataAll(:,:,ii);
    % Add zeros to account for latency
    dummyData = [dummyData;zeros(f.Latency,numChans)];
    y = squeeze(f(dummyData));
    % Remove latency segment
    y = y(f.Latency+1:end,:);
    filtAudio2(:,:,ii) = y;
    reset(f)
end
t2 = toc % 1.6368

% Compare results.
f1 = filtAudio;
f2 = filtAudio2;
norm(f1(:)-f2(:))

[jibrahim]

to your question on latency, yes, it would matter ina real-time context, where you are passing short frames of audio tot he object, and where a long latency might not be acceptable. in your context, since you process the entire audio file in one shot and then postprocess the output, it does not matter. you can deal with it by zero padding and then removing the unwanted signal part, as highlighted in my latest code snippet.

[Paul]

@jibrahim

Can you give any insight into why dsp.FrequencyDomainFIRFilter is so much faster than filter for these examples? My experience is that filter is blazingly fast, though I've never had the need to use it with a 240000 point impulse response.

[jibrahim]

Hi Paul,

dsp.FrequencyDomainFilter performs filtering in the frequency domain using fft-based overlap-add. As the filter length increases, fft becomes much faster than perfoming the filtering in the time domain, which is what the filter function does. The filters in this code are a few seconds long, so hundreds of thousands of samples, which explains the great difference in speed.

[Michael]

Thanks for the added explanation @jibrahim. @Paul I slipped up when I added dummy data into the example code - the impulse responses in question are not actually that long!