# A Narrow Bandpass Filter in Octave or Matlab

The design of a very narrow bandpass FIR filter, coded in either Octave or Matlab, can prove challenging if a computationally-efficient  filter is required. This is especially true if the sampling rate is high relative to the filter's center frequency. The most obvious filter design methods, using either window-based or Remez ( Parks-McClellan ) functions, can easily result in filters with many thousands of taps.

The filter to be described reduces the computational effort (and thus processing time) by a factor greater than 20 for Matlab, or 100 for Octave, compared to a single stage filter of similar performance.  For instance, at a 48 kHz sample rate, passbands from 1 Hz to over 2000 Hz, with excellent shape factor are possible with minimal design effort.

Fig 1.    The filter has complex down-conversion to baseband followed by filtering then up-conversion.  See Appendix A for a more detailed view of signal data within the filter.

THEORY OF OPERATION

Take a look at the filter block diagram above (Figure 1), which is taken from a filter design document by Momentum Data Systems[1].

This document gives the formula for the bandpass filter of Figure 1 as follows :

y0(ω) = (h(ω - ω0) + h(ω + ω0)) x0(ω)

The input signal, x(n), is processed as a vector of samples (eg from a '.wav' file), which are complex down-converted to baseband by two multipliers. Sine and cosine waveforms at the filter center frequency have the same sample length as the input signal vector. They are multiplied by the input to give baseband signals to each of the lowpass filters.

The two lowpass filters are identical, and use the MIFIR (Matrix IFIR) configuration which was described in my previous blog. The cut-off frequency of both filters is equal to half that of the desired bandpass width. After filtering, the complex up-conversion process converts the signal back to the input frequency.  For more interesting information on complex (quadrature) techniques, I can recommended this blog.
Note that the sample rate remains constant throughout the filter. In other words, there is no downsampling and upsampling.
There are just five parameters to be set which are as follows:

Bandpass width,  transition region width,  center frequency,  taps1,  taps2.

The taps1 parameter determines the stopband attenuation and passband ripple, while the taps2 value controls suppression of spurious products in the stopband.  The taps1 value will normally be greater than five times that of taps2.  It is suggested that the values given in the example (taps1 = 89 and taps2 = 15) are used as a starting point, and small changes made to tailor the response.

CODE WITH EXAMPLE PARAMETERS

The code listing below features parameters for a bandpass filter with a passband width of 10 Hz and ripple of less than 1 dB. The transition region width is 13 Hz wide with greater than 60 dB attenuation outside this range.
So what’s all this code doing ?   Well, despite the rather long software listing below, once the filter coefficients are calculated from the ‘remez’ routines, the filter process requires just fourteen lines.  Most of the remainder is to show plots of the results.
Note that there are a few combinations of input parameters for which the ‘remez’ algorithm may not converge.  This can be corrected by small changes to either the input parameters, or the value of parameters taps1 and/or taps2.
The software runs without modification in Octave version 5.1 or Matlab R2016b.
Figure 2 shows the filter response obtained with default values for taps1 and taps2.  View more detailed plots, including an FFT of the chirp output signal, by running the code.

Fig 2.  The filter response here uses the default design values of 48 kHz sample rate, 850 Hz center frequency, 10 Hz bandwidth, 13 Hz transition region, with values for taps1 and taps2 of 89 and 15 respectively.

The filter processing time is the time taken to process the input signal (eg the 5 second chirp signal in the example), between locations (A) and (H) of Figure 1.
This shows up in the Octave or Matlab command window, and at the top of each plot.
If the 'comparison_filter' variable near the top of the listing is set to 'true', then the filter processing time can be compared with that of a conventional single stage filter. The result is displayed in the Octave or Matlab command window.  Figure 3 shows the frequency response plot for the 50,000 tap single stage filter.

Fig 3.  The single stage 50000 tap comparison filter has a much longer processing time – compare this with Figure 2 (both run in Octave on Windows 10).

CONCLUSION

The efficient bandpass filter described has been designed for easy implementation in either Octave or Matlab. Comments and suggestions for adapting and/or improving the code are most welcome.

%
%                       FIR narrow bandpass filter
%                  ==========================

%                           Paul Lovell   June 2020

%   This program runs in Octave 5.1 using the 'remez' function
%   (from the Octave-Forge 'signal' package), or in Matlab.
%
%   Specification for this bandpass filter example :
%   Center frequency:               850 Hz
%   Passband:                            10 Hz at < 1dB
%   Transition Region:              13 Hz at > 60dB
%   Sample rate:                         48 kHz
%
%   The taps in the first and second stages of the two lowpass
%   filters can be set to any number > 3, as required.
%   Higher stopband attenuation, and/or a sharper transition
%   band will require more taps.
%

%  ***********************  Input Parameters  **********************

%   Set the bandpass filter input parameters as required
%
clear all;
Fs = 48000;            % Set sample rate in Hz
passband = 10;         % Set passband width Hz
stopband = 13;         % Set stopband limit Hz
cfr = 850;             % Set the center frequency
taps1 = 89;            % Set taps for first stage MIFIR filter
taps2 = 15;            % Set taps for second stage MIFIR filter

%  Set the test signal duration as required.
signal_block_drn = 5;    % test signal duration in seconds

comparison_filter = false;
%  This is used to check the relative processing duration and
%  compare it with a window method single stage filter

dt = 1/Fs;
tw = 0:dt:signal_block_drn;

psb = passband/2;
stb = stopband/2;

%  Test Chirp Signal input to the filter
%  Replace this with your own input signal as required

chfmin = 0; chfmax = Fs/2;
sig = chirp(tw,chfmin,signal_block_drn,chfmax)';
siglen = length(sig);

%  *****************  Down/up-conversion oscillators  *************

sinosc = sqrt(2) * sin(2*pi*cfr*tw)';
cososc = sqrt(2) * cos(2*pi*cfr*tw)';

%  *******  Calculate filter passband and stopband parameters  *******

m2 = round(Fs/(6 * (psb+stb)));      % interpolation for filter 2
m1 = 3 * m2;                         % interpolation for filter 1

pf = double(2 * psb * m1)/Fs;
sf = 1-pf;

f = [0 pf sf 1];                     % MIFIR-1 design parameters
b1 = remez(taps1-1, f, [1 1 0 0], [1 1]);

f = [0 0.199 0.451 1];               % MIFIR-2 design parameters
b2 = remez(taps2-1, f, [1 1 0 0], [1 8]);

scale = 1.02/sqrt(m1^2*m2);          % scaling factor
b1 = b1*scale; b2 = b2*scale;        % scale coefficients

bsig = zeros(length(sig),2);

%  *******************  Signal down-conversion  ****************

tstart = tic;

bsig(:,1) = sig .* cososc;
bsig(:,2) = sig .* sinosc;

% *********  Five stage baseband lowpass filters ********
%
for ix = 1:2

sig2 = [bsig(:,ix); zeros(mod(-siglen,m1),1)]; % zero padding
sg1 = filter(b1,1,reshape(sig2,m1,[])')';      % 1st stage (MIFIR)
sg2 = filter(b2,1,reshape(sg1,m2,[])')';       % 2nd stage (MIFIR)

sg3 = reshape(sg2,[],1);

sg3 = cumsum(sg3-[zeros(m1,1); sg3(1:end-m1)]); % 3rd stage RRS
sg3 = cumsum(sg3-[zeros(m1,1); sg3(1:end-m1)]); % 4th stage RRS

sg4(:,ix) =  cumsum(sg3-[zeros(m2,1);sg3(1:end-m2)]);
% 5th stage RRS
end

%  *******************  Signal up-conversion  ****************

sg5 = sg4(1:siglen,1) .* cososc;
sg6 = sg4(1:siglen,2) .* sinosc;
sig_out = sg5 + sg6;

timefilt =  num2str(toc(tstart));
disp(['Processing with dual 5 stage filter  ' ...
timefilt '  seconds']);

%  ******** FFT routines for frequency response plots **********

FFTsize = 2^(17+floor(log10(Fs/psb))); % Adjust FFT to passband
fp = Fs/1000;
b1 = b1/scale; b2 = b2/scale;         % Rescale to display MIFIRs
b3 = ones(m1,1)*scale; b4 = b3;       % Rescale to display RRS
b5 = ones(m2,1);
[H1,F1] = freqz(b1,1,FFTsize,fp);     % Spectrum of prototype filter
b1u = upsample(b1,m1);
[H1u,F1] = freqz(b1u,1,FFTsize,fp);   % Spectrum of matrix filter 1

[H2,F2] = freqz(b2,1,FFTsize,fp);     % Spectrum of prototype filter
b2u = upsample(b2,m2);
[H2u,F2] = freqz(b2u,1,FFTsize,fp);    % Spectrum of matrix filter 2

[H3,F3] = freqz(b3,1,FFTsize,fp);     % Spectrum of filter 3 (RRS)
[H4,F4] = freqz(b4,1,FFTsize,fp);     % Spectrum of filter 4 (RRS)
[H5,F5] = freqz(b5,1,FFTsize,fp);     % Spectrum of filter 5 (RRS)
[H6,F6] = freqz(sig_out,1,FFTsize,fp);

H = H1u.*H2u.*H3.*H4.*H5; % The overall filter frequency response

%  ********************  Graphic Display Function ******************

% amplitude vs frequency plots to show passband and stopband
% x-axis may be changed with a different plot width parameter
% this is the first of the four parameters for the plots below

filter_plot(1,:) =  [2 10 10 -80] ;
filter_plot(2,:) =  [20 10 10 -80] ;
filter_plot(3,:) =  [2 0.5 0.2 -2] ;
filter_plot(4,:) =  [85 10 10 -80] ;

H = vertcat(flipud(abs(H)),abs(H));    % Format the data
H = H(1:2:length(H));

for fig_num = 1:4

figure(fig_num);

plot(20*log10(H),'-r', 'linewidth', 2);
% Calculate axes based on FFT size
pg = psb * filter_plot(fig_num,1);
dmax = (FFTsize/2) + ((FFTsize * pg)/Fs);
dmin = (FFTsize/2) - ((FFTsize * pg)/Fs);
if cfr<pg
dmin = dmin + ((dmax-dmin)*(pg-(cfr/2))/(pg*2));
cp = 0; cr = Fs/2; cq = cr/20;
else
cp = cfr-pg; cr = cfr+pg; cq = pg/10;
end;
amax = filter_plot(fig_num,2);
intv = filter_plot(fig_num,3);
amin = filter_plot(fig_num,4);
axis([dmin,dmax,amin,amax]); grid on;
xlabel('Frequency, Hz'); ylabel('Amplitude, dB');

set (gca, 'xtick',[dmin:(dmax-dmin)/20:dmax]);
set (gca, 'xticklabel',[cp:cq:cr]);
set (gca, 'ytick',[amin:intv:amax]);
set (gca, 'yticklabel', [amin:intv:amax]);

title_str =(['Overall Frequency Response  -  ', ...
'Filter processing time  =  ',timefilt,' seconds']);
title(title_str);

end

figure(5);               %  FFT plot of filter output signal

plot(20*log10(abs(H6)/sqrt(Fs*2)),'-r', 'linewidth', 2);

dmax = 1275 * FFTsize * 2/Fs;        % Calculate axes based on FFT size
dmin = dmax/3; amax = 10; intv = 10; amin = -80;
smax = 1275; smin = 425;
axis([dmin,dmax,amin,amax]); grid on;
xlabel('Frequency, Hz'); ylabel('Amplitude, dB');

set (gca, 'xtick',[dmin:(dmax-dmin)/20:dmax]);
set (gca, 'xticklabel',[smin:(smax-smin)/20:smax]);
set (gca, 'ytick',[amin:intv:amax]);
set (gca, 'yticklabel', [amin:intv:amax]);

title_str =(['FFT of filter output sequence -  ', ...
'Filter processing time  =  ',timefilt,' seconds']);
title(title_str);

if comparison_filter
% Compares response with single stage filter (Hamming window)
% for default parameters.

b_ref = fir1(50000,[0.035208 0.035625]);
[Href,Fref] = freqz(b_ref,1,FFTsize,Fs); % Frequency response
figure(6); plot(Fref,20*log10(abs(Href)),'-k');
axis([840,860,-80,10]); grid on; xlabel('Frequency, Hz');
ylabel('Amplitude, dB');
set (gca, 'xtick',[840:1:860]);
set (gca, 'xticklabel',[840:1:860]);

ref = tic;
ref_filter_out = filter(b_ref,1,sig);
timeref = num2str(toc(ref));
t_str = ([' Single 50000 tap filter  -  ', ...
'Filter processing time  =  ',timeref,'  seconds']);
disp([char(10),'Processing with ',t_str,char(10)]);
title(strcat('Response with ',t_str));

end

%  ******************************************************************
%
%  This software is provided "AS IS", without warranty. The author
%  shall not be liable for any claim, damages or other liability
%  arising from, or in connection with this software.

APPENDIX  A

Fig 4.  A detailed analysis of signals at different points in the filter. Complex down-conversion is followed by filtering and complex up-conversion.

REFERENCES  AND  ACKNOWLEDGEMENTS

My thanks to Richard (Rick) Lyons for assistance with providing the data in Appendix A, and also with Matlab testing.

[1] Momentum Data Systems, http://www.mds.com/wp-content/uploads/2014/05/multirate_article.pdf, Page 8.

[ - ]
Comment by June 3, 2020

Very nice approach, especially because the coefficients for the baseband multiplications could be obtained from a numerically-controlled oscillator.

I could definitely see myself using this in the future!

[ - ]
Comment by June 3, 2020

Thanks therationalpi. There are plenty of possible hacks to the code to suit particular applications. Experiment and enjoy !

[ - ]
Comment by June 4, 2020

Hi Woodpecker,

Thanks for the article. So in effect you are saying: bring the band to dc, filter it then put it back where it was in frequency centre.

While this solution may be considered for software based dsp but not so for ASIC/FPGA as you increase design modules and effort to two mixers plus two filters.

By the way you stated that "The filter to be described reduces the computational effort (and thus processing time) by a factor greater than 20 for Matlab, or 100 for Octave, compared to a single stage filter of similar performance."

I hope Mathworks didn't read that. My understanding is that Octave is struggling to catch up with Matlab.

[ - ]
Comment by June 4, 2020

Hi kaz,

You're right, complex mix to baseband - then complex mix back up again.

I'm a little surprised you say that this would not be a good solution for ASIC/FPGA. However, one area you would have to watch is memory usage.

By the way, Mathworks doesn't have to worry. What I found is that the reference filter runs 5 times faster in Matlab than in Octave. And the new method described runs considerably faster than both of them :-)

[ - ]
Comment by April 30, 2021

This is a great article and very useful code. I am loathe to use it in projects though without a known license. Would you please consider releasing the code in this article under MPL or other permissive license?

[ - ]
Comment by May 1, 2021

Hi Nodal,

The code in this article should be considered experimental, and for personal use only.  It is not intended to be used in a commercial application.

Note that I have published an update, which has a number of improvements.

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.