Hi! I am going through a MATLAB based textbook for time frequency analysis of Neural data (MEG/EEG time series data) and got to a chapter which goes through filter construction. It’s primarily going over the usage of the functions firls (Least-squares linear-phase FIR filter design) and fir1 (Window-based FIR filter design) (see: docs and more docs). As a newcomer to this, I couldn’t find similar functions and was wondering if I could be pointed into the right direction. I have seen the filter docs in DSP.jl but they are a little difficult for me to understand. So
Which functions should I be using in the Julia ecosystem to do the same thing.
If the approach is different: how do they differ from Least-squares linear-phase FIR filter design and Window-based FIR filter design.
I should add I am looking for this because I’m doing the material in Julia, not in MATLAB. So I’m interested in learning and implementing, not comparing the methods between languages or anything.
(Sorry for the direct tag @ararslan, but I was informed you might be able to help with this, and I wanted to make a thread rather than a direct message.)
DSP.jl provides several ways to design filters. In your application, I suspect that the specific filtter design method is not very important. What you need to know is: the filter type (low pass, band pass, high pass); the sampling frequency; the filter order; and the cutoff frequency.
To design a filter, DSP.jl requires the specification of a design method, in addition to the parameters given above. For example, let’s say you want a low-pass filter with sampling frequency 200, cutoff frequency 50, and order 21. For low pass filters, Butterworth is a good, easy design method. Then, you can run:
using DSP
fs = 200; # sampling frequency
fc = 50; # cutoff frequency
order = 21; # filter order
filter = digitalfilter(Lowpass(fc, fs=fs), Butterworth(order));
This produces a filter with the following magnitude response:
Now let’s say you want to filter a signal x. All you need to do is run
You are correct, of course: indeed, the phase is slightly non-linear (the red line is the cutoff frequency):
Unfortunately DSP.jl does not support Bessel filters AFAICT. If perfectly linear phase is required, probably the best option is to use remez. Here’s a similar low-pass filter. Note that a larger order is required, and it has bandpass ripple (as often in filter design, it’s a matter of tradeoffs):
Thank you for your replies! They really helped me getting started. One of the things that author does is define a specific shape of the filter (I think to make a smoother transition?) with the following Matlab code:
So ffrequencies defines the frequencies of interest, normalized to the nyquist frequency. And then the shape of the filter is defined by idealresponse. Is there a way to make a smoother transition (assuming this is the reason) like this in Julia?
Also the filter order here comes out to something like 192, this doesn’t work with Butterworth(192). What exactly is this value? And is it the same as the filter order for butterworth and the “order- n FIR filter” (sorry if this latter part may be beyond the scope of this discussion)
@ElectronicTeaCup That code can be adapted to work with DSP.jl’s remez, producing a very similar filter. See my example above, and follow up if you need more help.
Thank you for the code and information. I checked out the documentation, and did the following to reproduce something similar to the MATLAB code I was using.
center_freq = 20 # in Hz
filter_frequency_spread = 6 # Hz +/- the center frequency
transition_width = 0.2
ffrequencies = [
0,
(1-transition_width)*(center_freq-filter_frequency_spread),
(center_freq-filter_frequency_spread),
(center_freq+filter_frequency_spread),
(1+transition_width)*(center_freq+filter_frequency_spread),
nyquist,
] /nyquist
idealresponse = [ 0, 0, 1, 1, 0, 0, ]
# Order hard set to 200 as in MATLAB code
filterweights = remez(
200,
[
(ffrequencies[1], ffrequencies[2]) => 0,
(ffrequencies[3], ffrequencies[4]) => 1,
(ffrequencies[5], ffrequencies[6]) => 0
],
Hz = EEG["srate"]/nyquist
)
And then to take a look at what’s going on with a fft
filterweights = (firls(200,ffrequencies,idealresponse));
# compute its power spectrum
fft_filtkern = abs.(fft(filterweights))
fft_filtkern = fft_filtkern./maximum(fft_filtkern) # normalized to 1.0 for visual comparison ease
hz_filtkern = range(0,nyquist, length=101); # list of frequencies in Hz corresponding to filter kernel
The problem you’re seeing is that the filter you want cannot be created within the parameters you’ve specified. In broad terms:
Narrow filters (relative to the Nyquist frequency), require large orders and/or wide transition bands.
Very small transition bands also require large filter orders.
Bandpass filters with lower cutoff frequency close to 0, or higher cutoff frequency close to Nyquist, are in general also hard to design.
The degrees of freedom you have available are: filter bandwidth, filter order, and transistion bandwidth. Relaxing one can improve the others. You may need to try a few combinations before you get a filter you like.
I gave it a shot with the bandwidth you need (20 Hz). I relaxed the transition bandwidth to 2 Hz and increased the order to 300. I also asked remez to work harder by increasing maxiter. This is what I got:
using DSP
fs = 250.0 # sampling frequency
fn = fs/2.0 # Nyquist frequency
B = 20.0 # filter bandwidth
fc = 20.0 # center frequency
fl = fc - B/2. # lower cutoff
fh = fc + B/2. # higher cutoff
tr = 2.0 # transition bandwidth
order = 400 # filter order
# calculate coefficients
coeff = remez(order, [(0, fl-tr) => 0,
(fl, fh) => 1,
(fh+tr, fn) => 0],
Hz = fs,
maxiter = 50)
# find and plot magnitude response
f = range(0, fn, length = 250)
fz = freqz(PolynomialRatio(coeff,[1.0]), f, fs)
mag = abs.(fz)
plot(f, mag)
Ah right, it’s majorly the asymetric transition zones in my example that seems to be causing the issue. Thank you so much, I didn’t have a clue that these parameters would make all the difference!
I do not know why exactly the transition zones differ on each end in the text I’m following, but I want to try to get the same outputs as the rest of the analysis. Do you think it would be possible to get a transition like the transition regions of the first block of code?
In my opinion, you’re hitting the limits of what remez can do. Unfortunately, firls and other more modern design algorithms don’t exist in DSP.jl yet, and FIRWindow does not allow specifying the transition bands.
Note that Matlab also struggles with this filter. Using firls:
fir1 does not support specifying transition bands (this command is equivalent to DSP.jl’s FIRWindow design method). The best results are obtained with fir2, but note that the transition bands seem quite a bit wider than specified:
What I would do in your case is to use remez to design a filter with narrow transition bands, flat passband, and good (40 dB or more) rejection, and move on If you absolutely want to use the same filter as the textbook, you can generate the filter coefficients in Matlab and copy/paste them to your Julia program. If you don’t have access to Matlab, send me the code and I’ll send you the coefficients.
I hope I can find the time soon to implement firls and other design algorithms for DSP.jl, if nobody else beats me to it…
I think I’ll follow your advice, and move on working with what I’ve got—thanks for the MATLAB offer btw. I really don’t want to go down the using MATLAB with Julia path—since the whole point was to learn Julia in this context. Thanks for all the support! I hope you end up porting firls
If you absolutely want to use the same filter as the textbook, you can generate the filter coefficients in Matlab and copy/paste them to your Julia program. If you don’t have access to Matlab, send me the code and I’ll send you the coefficients.
