Update to a Narrow Bandpass Filter in Octave or Matlab
Following my earlier blog post (June 2020) featuring a Narrow Bandpass Filter, I’ve had some useful feedback and suggestions. This has inspired me to come up with an updated version, incorporating the following changes compared to the earlier one :
- Simpler code in Octave or Matlab
- Float32 precision replaces float64
- Faster processing by a factor of at least 4 times
- Easier setup of input parameters
- Normalized signal output level
A new experimental version in Python and Tensorflow running on Google Colab is also available here, more details at the end of this blog post.
OPERATION AND UPDATE DETAILS
For new readers, I suggest reading my previous blog which explained the theory of operation. As a reminder, Figure 1 shows the basic configuration with complex down-conversion followed by baseband lowpass filtering and complex up-conversion.
Fig 1. The basic filter with input signal x(n) and output y(n). [1] Note that the sample rate throughout the filter remains constant.
Each of the two baseband filters of Figure 1 now has a single matrix IFIR filter followed by three identical RRS (moving sum) filters, each of length m1 (the IFIR expansion factor). The m1 expansion value is calculated automatically in the software. Note that the two baseband filters are identical.
Fig 2. Detail of a single lowpass baseband filter from Figure 1.
Filter design parameters include the number of taps and passband weight which are applied to the ‘remez’ filter design function (‘firpm’ may be used in Matlab). These function as for a conventional single stage filter, in that more taps give less passband ripple and higher stopband attenuation.The ‘weight_sb’ input parameter determines the ratio of stopband attenuation to passband ripple. Experimenting with both these values will give a useful feel for the capabilities of the filter.
The IFIR filter expansion value ‘m1’ is set by the software based on the sample frequency, number of taps, and the filter bandwidth.
The frequency response is illustrated using a chirp signal as the filter input. To obtain an accurate plot requires a minimum ‘signal block duration’ input parameter determined by the number of taps, and the passband. Try using: sbd = taps/(passband*2) which gives 5 seconds for the example parameters of 10 Hz passband and 100 taps.
An external signal can replace the chirp test signal. This can have any of the following data formats:16, 24 or 32 bit unsigned integer, float32 or float64. Input data is converted to float32 for processing. The code has been tested with ‘.WAV’ files for input and output.
CODE WITH EXAMPLE PARAMETERS
Once again, the code listing below includes example parameters for a bandpass filter with a passband of 10 Hz and ripple less than 1 dB. At a bandwidth of 13 Hz or greater the stopband attenuation exceeds 60 dB. Testing has been successful using a very wide range of sampling rates, bandwidths, and passband/stopband specifications.
The results in Figure 3 shows FFT plots generated using a filter with the example parameters. Running the code will also give a zoomed-in plot of the passband response.
Fig 3. Plots show filter frequency response after running the program in Octave with the example parameters from the code below.
CONCLUSION
The efficient bandpass filter described has been designed for easy implementation, and runs without modification in GNU Octave or Matlab R2019b.
A new experimental version running Python and Tensorflow in Google Colab is also available here, which can be run by those of you with a Google account. To enable a GPU in Colab which considerably increases the processing speed :
Select :
- Menu> Runtime>Change runtime type.
- Change hardware acceleration to GPU.
- Click save.
I’m fairly sure the Tensorflow code could be made more efficient, in particular to make better use of the GPU.
Comments and suggestions for adapting and/or improving the code of either version are most welcome.
% % FIR narrow bandpass filter ver 2.0 % =================================== % % Paul Lovell March 2021 % This program runs in Octave 6.1 using the 'remez' function % (from the Octave-Forge 'signal' package), or in Matlab. % It is an updated version of the design first published on % dsprelated.com in June 2020. % % 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 stage lowpass filter can be set % to any number > 30, as required. % Higher stopband attenuation, and/or a sharper transition % band will require more taps. % The weight of the stopband can also be increased or % decreased relative to the passband. This results in % changes to passband ripple and stopband attenuation. % % *********************** 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 taps = 100; % Set taps for first stage matrix IFIR filter sb_weight = 20; % Weight of stopband attenuation relative % to passband ripple (passband weight = 1) % Set the test signal duration in seconds as required. % The minimum value here is approx : taps/(passband*2) % (eg 5 seconds for a 10Hz passband filter with 100 taps) signal_block_drn = 5; % test signal duration in seconds profile_on = false; % The profile option works in Octave but may need to be % changed for Matlab 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 = cfr * 2; dt = 1/Fs; tw = (0:dt:signal_block_drn); sig = chirp(tw,chfmin,signal_block_drn,chfmax); siglen = length(sig); %% % Enable this section for input from single channel .WAV file %% %% [sig, Fs] = audioread("c:/dsp/input.wav"); % alternative input %% sig = sig'; siglen = length(sig); %% dt = 1/Fs; signal_block_drn = siglen/Fs; %% tw = (0:dt:signal_block_drn); % ***************** Down/up-conversion oscillators ************* sinosc = single(sqrt(2) * sin(2*pi*cfr*tw)); cososc = single(sqrt(2) * cos(2*pi*cfr*tw)); % ******* Calculate filter passband and stopband parameters ******* m1 = round(Fs/(psb * round(9+(taps/30)))); % interpolation for % first stage pf = 2 * psb * m1/Fs; sf = pf * stopband/passband; f = [0 pf sf 1]; % MIFIR design parameters b1 = remez(taps-1, f, [1 1 0 0], [1 sb_weight]); % if size(b1,1)>1 b1=b1'; end; b1 = single(b1/(m1^3)); % scale coefficients profile clear; if profile_on profile on; end; tstart = tic; % ********* 2 x Four stage baseband lowpass filters ******** zk = zeros(1,mod(-siglen,m1)); zr = zeros(1,m1); zt = siglen-m1; sg1 = reshape([sig.*cososc, zk],m1,[]); % Signal down-conversion sg2 = reshape(conv2(sg1,b1),1,[]); sg2 = sg2(1:siglen); % 2D convolution sg3 = cumsum(sg2-[zr, sg2(1:zt)]); % 2nd stage RRS sg4 = cumsum(sg3-[zr, sg3(1:zt)]); % 3rd stage RRS sg5a = cumsum(sg4-[zr, sg4(1:zt)]); % 4th stage RRS sg1 = reshape([sig.*sinosc, zk],m1,[]); % Signal down-conversion sg2 = reshape(conv2(sg1,b1),1,[]); % 2D convolution sg2 = sg2(1:siglen); sg3 = cumsum(sg2-[zr, sg2(1:zt)]); % 2nd stage RRS sg4 = cumsum(sg3-[zr, sg3(1:zt)]); % 3rd stage RRS sg5b = cumsum(sg4-[zr, sg4(1:zt)]); % 4th stage RRS % ******************* Signal up-conversion **************** sig_out = (sg5a .* cososc)' + (sg5b .* sinosc)'; tf = toc(tstart); profile off; if profile_on profshow(); end; timefilt = num2str(tf); disp(['Processing with dual 5 stage filter ' ... timefilt ' seconds']); ## % Enable this section for output to single channel .WAV file ## audiowrite("c:/dsp/output.wav",sig_out',Fs); ## return % ******************** Graphic Display Function ****************** % Amplitude vs frequency plots from FFT of output signal. FFTsize = 2^(16+floor(log10(Fs/psb))); % Adjust FFT to passband fp = Fs/1000; pb = passband; [H,F] = freqz(sig_out,1,FFTsize,fp); % x-axis limits may be changed using the first two % of the five parameters for the plots below (Hz) filter_plot(1,:) = [cfr-pb cfr+pb 10 10 -80] ; filter_plot(2,:) = [cfr-pb cfr+pb 1 0.2 -2]; filter_plot(3,:) = [0 cfr*2 10 10 -80] ; for fig_num = 1:3 figure(fig_num); cp = filter_plot(fig_num,1); cr = filter_plot(fig_num,2); cq = (cr-cp)/10; plot(20*log10(abs(H/max(H))),'-r', 'linewidth', 2); fg = 'Response of filter output signal (FFT) - '; dmax = cr * FFTsize * 2/Fs; % Calculate axes based on FFT size dmin = cp * FFTsize * 2/Fs; amax = filter_plot(fig_num,3); intv = filter_plot(fig_num,4); amin = filter_plot(fig_num,5); axis([dmin,dmax,amin,amax]); grid on; xlabel('Frequency, Hz'); ylabel('Amplitude, dB'); set (gca, 'xtick',[dmin:(dmax-dmin)/10:dmax]); set (gca, 'xticklabel',[cp:cq:cr]); set (gca, 'ytick',[amin:intv:amax]); set (gca, 'yticklabel', [amin:intv:amax]); title([fg,'Filter processing time = ',timefilt,' seconds']); 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.
REFERENCES AND ACKNOWLEDGEMENTS
My thanks to Richard (Rick) Lyons for assistance with editing and Matlab testing.
[1] Momentum Data Systems, http://www.mds.com/wp-content/uploads/2014/05/multirate_article.pdf, Page 8.
- Comments
- Write a Comment Select to add a comment
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.
Please login (on the right) if you already have an account on this platform.
Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: