Quickest way to do a bandpass filter?

I have a huge dataset of .wav files that I would like to bandpass filter.

Is the following the most time efficient way to do it..?

This filter takes about 2 minutes to run so it's quite time consuming with such a huge dataset.

[xbit, fs]=audioread(filename, [1*fs,120*fs]); %read in file from 1-120 secs

xbit=detrend(xbit);%removes DC offset

tic
xbit_filt=bandpass(xbit,[50 24000],fs);
toc

Answer 1

[Mathieu NOE]

Hello Louise

welcome back and happy new year !

I made some speed comparison with your code and two other options using a butterworth digital filter , then using the function filter or filtfilt (filtering twice , forward and backward in time, for phase distorsion cancellation)

see my suggestion are much faster than using the bandpass function as in your code

I wonder why you put the low pass frequency at 24 kHz as this is already above the human hearing upper frequency limit...

at what sampling frequency are your wav files ?

with the my small audio file sampled at 11025 Hz , the differences are significant :

• your code : Elapsed time is 1.574109 seconds.
• my code with filter : Elapsed time is 0.016102 seconds.
• my code with filtfilt : Elapsed time is 0.042606 seconds.

filename = 'test_voice_mono.wav'  ;
% NB this file is sampled at 11025 Hz, so the bandpass filter cut off
% frequencies are modified  (LPF : 24kHz => 2.4kHz)

[xbit, fs]=audioread(filename); %we don't know yet what if fs so I cannot define range = [1*fs,120*fs] as in your original code
samples = length(xbit);
start = 1*fs;
stop = min(120*fs,samples);
xbit=xbit(start:stop); % extract data for t = 1-120 secs

xbit=detrend(xbit);%removes DC offset

% your code
tic
xbit_filt=bandpass(xbit,[50 2400],fs);
toc


% good old butter bandpass filter (applied once)
tic
N = 2; 
[B,A] = butter(N,[50 2400]/(fs/2));
xbit_filt = filter(B,A,xbit); % NB filter applies the filter once ! 
toc


% good old butter bandpass filter (applied twice)
tic
N = 2; 
[B,A] = butter(N,[50 2400]/(fs/2));
xbit_filt = filtfilt(B,A,xbit); % NB filtfilt applies the filter twice ! 
toc

[Louise Wilson]

Thank you Mathieu! Happy New Year :-)

Sorry, I forgot the sampling rate - 144,000 kHz.

Do you mean upper limit with regard to 24kHz? I am looking at sound in the ocean so not so much focused on human hearing.

Excited to see there is a possibility to speed up! I'll try this now. What is the motivation behind applying the filter twice?

Thanks so much for your help.

[Mathieu NOE]

Now you can build a for loop to apply this code to all wav files in your directory

If for some reasons you need the files to be sorted (really) , see this FEX submission

Natural-Order Filename Sort - File Exchange - MATLAB Central (mathworks.com)

otherwise you can skip / comment that line : S = natsortfiles(S);

%% define path
fileDir = pwd; % or your specific path / folder

S = dir(fullfile(fileDir,'*.wav')); % get list of data files in directory  
S = natsortfiles(S); % sort file names into natural order (what matlab dir does not well) , see FEX : 
           %(https://fr.mathworks.com/matlabcentral/fileexchange/47434-natural-order-filename-sort)

               
% good old butter bandpass filter (applied once)
tic

               
%% Loop inside folder
for k = 1:length(S) % read data in specified sheet
    fileName = S(k).name;
    [xbit,fs] = audioread(fullfile(fileDir, fileName)); % or use a structure (S(k).data ) to store the full data structure
    dt = 1/fs;
    [samples,channels] = size(xbit);

    start = 1*fs;
    stop = min(120*fs,samples);
    xbit=xbit(start:stop); % extract data for t = 1-120 secs
    time = (start:stop)*dt; % time vector 

    xbit=detrend(xbit);%removes DC offset

    % good old butter bandpass filter (applied once)
    N = 2; 
    [B,A] = butter(N,[50 2400]/(fs/2));
    xbit_filt = filter(B,A,xbit); % NB filter applies the filter once ! 
    
    % add your own code here (save xbit_filt to wav file ?)
    

end
        
toc

[Louise Wilson]

Sorry Mathieu, I see you explained the reason to filter twice :-)

[Mathieu NOE]

Louise

Differences between filter and filtfilt

filter applies the defined filter (B,A coefficients) on the data only once (forward in time).

filtfilt applies the defined filter (B,A coefficients) on the data twice (forward and backward in time). this way you have no phase distorsion (filter lag) . This is off course only doable "off line" at it is not causal !

help filtfilt

filtfilt Zero-phase forward and reverse digital IIR filtering.

After filtering in the forward direction, the filtered sequence is then

reversed and run back through the filter; Y is the time reverse of the

output of the second filtering operation. The result has precisely

zero phase distortion, and magnitude modified by the square of the

filter's magnitude response.

If your target is to limit the signal to 24 kHz , why not simply decimate the signal with decimation factor 3 ?

this way you remove what is above 24 kHz and you reduce also the size of the data by same factor 3

function odata = decimate(idata,r,nfilt,option)
%DECIMATE Resample data at a lower rate after lowpass filtering.
%   Y = DECIMATE(X,R) resamples the sequence in vector X at 1/R times the
%   original sample rate.  The resulting resampled vector Y is R times
%   shorter, i.e., LENGTH(Y) = CEIL(LENGTH(X)/R). By default, DECIMATE
%   filters the data with an 8th order Chebyshev Type I lowpass filter with
%   cutoff frequency .8*(Fs/2)/R, before resampling.

see end of my code below

• your code : Elapsed time is 1.426317 seconds.
• my code with filter : Elapsed time is 0.002231 seconds.
• my code with filtfilt : Elapsed time is 0.002717 seconds.
• my code with decimate : Elapsed time is 0.008000 seconds.

filename = 'test_voice_mono.wav'  ;
% NB this file is sampled at 11025 Hz, so the bandpass filter cut off
% frequencies are modified  (LPF : 24kHz => 2.4kHz)

[xbit, fs]=audioread(filename); %read in file from 1-120 secs
[samples,channels] = size(xbit);
start = 1*fs;
stop = min(120*fs,samples);
xbit=xbit(start:stop); % extract data for t = 1-120 secs

xbit=detrend(xbit);%removes DC offset

% your code
tic
xbit_filt1=bandpass(xbit,[50 2400],fs);
toc

% good old butter bandpass filter (applied once)
tic
N = 2; 
[B,A] = butter(N,[50 2400]/(fs/2));
xbit_filt2 = filter(B,A,xbit); % NB filter applies the filter once ! 
toc


% good old butter bandpass filter (applied twice)
tic
N = 2; 
[B,A] = butter(N,[50 2400]/(fs/2));
xbit_filt3 = filtfilt(B,A,xbit); % NB filtfilt applies the filter twice ! 
toc

%% decimate by factor 3 
% NB : decim = 1 will do nothing (output = input)
tic
decim = 3;
if decim>1
    for ck = 1:channels
    xbitdecim(:,ck) = decimate(xbit(:,ck),decim);
    fs = fs/decim;
    end
end
samples = length(xbitdecim);
time = (0:samples-1)*1/fs;
toc

[Louise Wilson]

Thanks Mathieu for your brilliant and detailed answer.

I am just using the filter to get the 5th and 95th percentiles, median, and rms level of the file, then store those in an output. Got the loop sorted :-)

[Mathieu NOE]

ok

glad I could help !

all the best