Quickest way to do a bandpass filter?

Question

Louise Wilson 2023-1-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1890802-quickest-way-to-do-a-bandpass-filter

评论： Mathieu NOE 2023-1-9

采纳的回答： Mathieu NOE

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2023-1-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1890802-quickest-way-to-do-a-bandpass-filter#answer_1144797

在 MATLAB Online 中打开

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

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

Mathieu NOE 2023-1-9

在 MATLAB Online 中打开

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

Mathieu NOE 2023-1-9

在 MATLAB Online 中打开

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 2023-1-9

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 :-)