Multichannel signal filter with multichannel impulse response (efficiently)

Question

Michael 2024-10-23，10:08

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2161615-multichannel-signal-filter-with-multichannel-impulse-response-efficiently

评论： Michael 2024-10-24，8:13

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?

2 个评论
显示无隐藏无

Hitesh 2024-10-23，12:46

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 2024-10-23，14:15

在 MATLAB Online 中打开

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 2024-10-23，18:01

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2161615-multichannel-signal-filter-with-multichannel-impulse-response-efficiently#answer_1536420

在 MATLAB Online 中打开

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

jibrahim 2024-10-23，20:32

编辑：jibrahim 2024-10-23，20:39

在 MATLAB Online 中打开

(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 2024-10-24，1:59

编辑：jibrahim 2024-10-24，2:00

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 2024-10-24，8:13

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!

请先登录，再进行评论。

Multichannel signal filter with multichannel impulse response (efficiently)

2 个评论
显示无隐藏无

采纳的回答

6 个评论
显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Multichannel signal filter with multichannel impulse response (efficiently)

2 个评论 显示 无隐藏 无

采纳的回答

6 个评论 显示 4更早的评论隐藏 4更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

6 个评论
显示 4更早的评论隐藏 4更早的评论