0
$\begingroup$

I have a question regarding filtering (bandpass in my case) as a preprocessing step before a classification task.

The frequency band I am interested in lies between 0.5 - 22 Hz which results in a 6000th order bandpass filter due to the steep slope needed for the highpass. The phase response is linear and I use the filtfilt matlab function to get a zero phase output.

After filtering I try to classify on sample level, to determine when a certain state of two is present in my data.

My question now is whether the high filter order distorts the time sequence, since the output samples of the filter are a superposition of samples $n$ to $n-5999$?

$\endgroup$
3
  • $\begingroup$ 1 - What's your sampling frequency? 2 - Why do you feel you need to filter? 3 - Is it real-time or post processing? 4 - Have you considered using an IIR filter ? $\endgroup$
    – Ben
    Jan 3 at 18:56
  • $\begingroup$ Hi Ben. Thank you for your response! 1. My sampling frequency is 1kHz 2. I work with EEG signal which has several artifacts from head movements, eye blinks etc. in the higher frequency range. Also DC drifts are an issue. That's why I use a high pass filter as well. 3. The filtering is part of preprocessing before classification. At the moment real time is no issue. 4. Yes I have, but I want to avoid it because of stability issues. $\endgroup$
    – handsitter
    Jan 3 at 19:07
  • $\begingroup$ Since you use an FIR, you don't need to use filtfilt. You could simply remove the number of samples corresponding to the group delay of your FIR filter. $\endgroup$
    – Ben
    Jan 4 at 17:51
1
$\begingroup$

I would handle this as a chain of FIR filters with decimation to get much lower tap count FIR filters that don't need to be nearly as steep in the early stages.

Maybe using Blackman-Harris windowed filters for the first two stages (for superior alias fold-in rejection).

Then a Kaiser windowed filter for the 3 stage to get a flatter top and a sharper roll-off.

I guarantee you'll be using much less than 6000 taps total for the same filtering performance.

I've worked out fairly decent scheme in some Octave/MatLab code:

LPF -> decimate by 2 -> LPF -> decimate by 5 -> BPF

This gets you to a final 100 samples/second sample rate.

The three filters, that I came up with, have tap lengths of 9, 51, and 205 respectively.

The two low-pass anti-alias filters keep any alias fold-in into your band of interest at -40 dB or lower.

The final band-pass filter is -35 dB at DC, -6.3 dB at 0.5 Hz, -0.1 dB at 22.5 Hz, and -40 dB at 24.5 Hz.

enter image description here enter image description here enter image description here

Here's the Octave code.

Main script:

% First anti-alias filter, designed for decimation by 2 after the filter
Fs1 = 1000;  % [-500,500]
a = blackman_harris_lpf(1.0, Fs1, 100.0, 300.0, 55);
length(a)

% Second anti-alias filter, designed for decimation by 5 after the filter
Fs2 = Fs1/2;  % [-250, 250]
b = blackman_harris_lpf(1.0, Fs2, 50, 20.0, 45);
length(b)

% Band pass filter for the frequency band of interest
Fs3 = Fs2/5;  % [-50, 50]
c = kaiser_lpf2(1.0, Fs3, (23.5-0.5)/2, 2, 9);
radian_phase_inc = pi/(Fs3/2) * (24.5-0.5)/2;
c = 2*real(c .* exp(1j * radian_phase_inc .* [0:(length(c)-1)]'));
length(c)

% Get response curves for each filter for plotting
[ha,wa] = freqz(a, [1], 2*5*1000);
[hb,wb] = freqz(b, [1], 2*1000);
[hc,wc] = freqz(c, [1], 2*1000/5);

% Create plotting x-axis (aka frequency-axis) value arrays to
% fold freq axis to show aliasing into final passband.
% Scale to original Fs1
waturn = [wa(1:(2*5*1000/2))' wa((2*5*1000/2):-1:1)'] * 500/pi; % Fs = 1000, fold for decimation by 2
wbturn = [wb(1:(2*1000/5))'   wb((2*1000/5):-1:1)' ...
          wb(1:(2*1000/5))'   wb((2*1000/5):-1:1)' ...
          wb(1:(2*1000/5))'                         ] * 250/pi; % Fs = 500, fold for decimation by 5
wcturn = wc * 50/pi; % Fs = 100, no turns

clf;
plot( ...
     waturn, 20*log10(abs(ha)), ';1st filter, Blackman-harris LPF, 1000 sps folded to 500 sps;', ...
     wbturn, 20*log10(abs(hb)), ';2nd filter, Blackamn-harris LPF,  500 sps folded to 100 sps;', ...
     wcturn, 20*log10(abs(hc)), ';3rd filter, Kaiser BPF, 100 sps;'  ...
);
title('Filters Alias Folding Performance');
ylabel('Gain (dB)');
xlabel('Frequency (Hz)');
grid on;
hold off;

blackman_harris_lpf.m:

function b = blackman_harris_lpf(gain, Fs, Fc, W, atten_dB)
% B = blackman_harris_lpf(GAIN, FS, FC, W, ATTEN_DB)
%
% Design a Blackman-harris windowed low pass filter
%
%    GAIN Passband linear gain. E.g. 1.0
%
%    FS   Sampling frequency in Hz
%
%    FC   Cutoff frequency (the edge of the ideal LPF response before
%         before truncation and windowing of the sinc() pulse) in Hz.
%
%    W    Transition width in Hz
%
%    ATTEN_DB  stop band attenuation in dB

M = floor(atten_dB * Fs/(22.0*W)); % rule of thumb from fred harris
if mod(M,2) == 0  % ensure M is odd
   M = M + 1;
end

% Windowed, truncated sinc()
w = blackmanharris(M);
h = sinc([-(M-1)/2:1:(M-1)/2]*Fc/(Fs/2)).' .* w;

% Normalize passband response to 0 dB and apply gain
b = gain * h / sum(h);

kaiser_lpf2.m:

function b = kaiser_lpf2(gain, Fs, Fc, W, Beta)
% B = kaiser_lpf2(GAIN, FS, FT, W, BETA)
%
% Design a Kaiser windowed low pass filter
%
%    GAIN Passband linear gain. E.g. 1.0
%
%    FS   Sampling frequency in Hz
%
%    Fc   Cutoff frequency (the edge of the ideal LPF response before
%         before truncation and windowing of the sinc() pulse) in Hz.
%
%    W    Transition width in Hz
%
%    BETA Kaiser window parameter dictating the width of the main lobe
%         and the height of the stopband lobes.  Beta == 0 is equivalent
%         to a rectangular window.

M = floor((Beta/0.1102 + 8.7) * Fs/(22.0*W)); % rule of thumb from fred harris
if mod(M,2) == 0  % ensure M is odd
   M = M + 1;
end

% Windowed, truncated sinc()
w = kaiser(M, Beta);
h = sinc([-(M-1)/2:1:(M-1)/2]*Fc/(Fs/2)).' .* w;

% Normalize passband response to 0 dB and apply gain
b = gain * h / sum(h);
$\endgroup$
2
  • $\begingroup$ Hi Andy, thanks for answering. I will try your approach. It would be great, if I could end at 100 samples/second in the end, because downsampling is a step in my preprocessing pipeline anyway. Would it be ok to upsample before the filtering steps, so that I end at a specific sampling rate? $\endgroup$
    – handsitter
    Jan 4 at 10:43
  • $\begingroup$ @handsitter I edited my answer to provide an explicit solution for a 100 samples/second solution. $\endgroup$
    – Andy Walls
    Jan 10 at 0:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.