# High FIR order in time critical classification task

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$$?

• 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 ?
– Ben
Jan 3 at 18:56
• 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. Jan 3 at 19:07
• 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.
– Ben
Jan 4 at 17:51

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.

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

• 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? Jan 4 at 10:43
• @handsitter I edited my answer to provide an explicit solution for a 100 samples/second solution. Jan 10 at 0:29