Bandpass Filter (process_bandpass) - How does it work?

Vyte · March 13, 2025, 6:33pm

Hey everyone,

I wanted to make a post to ask about the Bandpass filter, specifically how its designed. I am running wPLI on some source EEG data so I used the Brainstorm GUI to generate a wPLI matrix and then I also wrote my own wPLI code on Matlab. The matrices are a perfect match (100% correlation), but I do not understand the filter and would like some explanation/advice.

So, to extract frequency within the desired frequency range (e.g. 1-4Hz Delta), I use the following Brainstorm function:

    DataAband = process_bandpass('Compute', my_data, Fs, (low_freq), (high_freq));

(my_data is regions x time, Fs is 250Hz sampling frequency, low/fast_freq are the band limits)

I looked here: https://neuroimage.usc.edu/brainstorm/Tutorials/ArtifactsFilter ...to find some explanation about this bandpass filter. I now understand this is a FIR type filter, and mostly likely 60dB stopband attenuation (set by default) but I know nothing else about it.

I tried to design this filter myself to avoid using the function as I am not sure about the specifics, but the results are not the same:

Fs = 250; % Sampling frequency (Hz)
low_cutoff = 1; % Lower cutoff frequency (1 Hz) % Define the bandpass filter cutoff frequencies
high_cutoff = 4; % Upper cutoff frequency (4 Hz)
nyquist = Fs / 2; % Normalize the cutoff frequencies by the Nyquist frequency (Fs/2)
low_cutoff_norm = low_cutoff / nyquist;
high_cutoff_norm = high_cutoff / nyquist;
order = 50; % Filter order (no idea what to set this to exactly)
b = fir1(order, [low_cutoff_norm high_cutoff_norm], 'bandpass'); % Filter order (again, no idea how to modify this)

I am wondering if there is a publication or something that explains all the details of this filter or if there is a way to know the specific parameters.
If I could understand how 'process_bandpass' works, then I wouldn't need to design the filter myself. I am just trying to build a good understanding of my own work and why this process works.

Any advice would be very welcome!

Healthy Regards,
Vyte

Raymundo.Cassani · March 13, 2025, 7:33pm

With such call, the function Compute in the process_bandpass.m file calls bst_bandpass_hfilter, where a Kaiser window FIR filter design is implemented

github.com/brainstorm-tools/brainstorm3

toolbox/process/functions/process_bandpass.m

c113225a8


      
          [x, FiltSpec, Messages] = bst_bandpass_hfilter(x, sfreq, HighPass, LowPass, isMirror, isRelax, [], TranBand, Method);

github.com/brainstorm-tools/brainstorm3

toolbox/math/bst_bandpass_hfilter.m

c113225a8


      
              fcuts(3:4) = [];
              mags(3) = [];
              devs(3) = [];
          end
          if (f_highpass == 0)  % User didn't want a highpass
              fcuts(1:2) = [];
              mags(1) = [];
              devs(1) = [];
          end
          
          % Generate FIR filter
          % Using Matlab's Signal Processing toolbox
          if bst_get('UseSigProcToolbox')
              [n,Wn,beta,ftype] = kaiserord(fcuts, mags, devs, 2);
              n = n + rem(n,2);  % ensure even order
              b = fir1(n, Wn, ftype, kaiser(n+1,beta), 'noscale');
              % Using Octave-based functions
          else
              [n,Wn,beta,ftype] = oc_kaiserord(fcuts, mags, devs, 2);
              n = n + rem(n,2);  % ensure even order
              b = oc_fir1(n, Wn, ftype, oc_kaiser(n+1,beta), 'noscale');

You can have all the details for the designed filter by clicking on the button View filter response in the GUI of the process Pre-process > Band-pass filter

Finally, you can use the Matlab function freqz and impz to plot the frequency and impulse responses of a digital filter. Remember that in FIR, the denominator a is always equal to 1.

b = fir1(500, [1, 4] ./ (sfreq ./ 2), 'bandpass'); 
[Hf,Freqs] = freqz(b, 1, 2^12, sfreq);
[Ht,t] = impz(b, 1, [], sfreq);
figure()
% Plot frequency response
subplot(2,1,1)
plot(Freqs, 20.*log10(abs(Hf)));
xlabel('frequency (Hz)')
xlabel('magnitude (db)')
% Plot impulse response
subplot(2,1,2)
plot(t, Ht);
xlabel('time (s)')
ylabel('amplitude')

If you are not sure what are the parameters for a FIR filter, you may want to check some fundamentals on digital signal processing. In brief, the larger the order, the sharper the filter, but the impulse response last longer. An the type of window define the filter lobes.