passive RC highpass magnitude FRF to MATLAB FIR coefficients

Hello,

I have a GRASS EEG LP511 amplifier. It has a high-pass filter that is described as passive with 12 dB/octave slope. I imagine it is something like two simple RC high-pass filters in series. I can measure the magnitude frequency response function easily and I know the cut-off frequency (-6dB at 3 Hz).

I would like to implement this as an FIR filter in MATLAB.

I am having trouble sorting out how to obtain the FIR filter coefficients to match my analogue high-pass filter.

Could anyone help? Please assume I don't have a strong enough background and will need thorough explanations. Many Thanks! DP

Answer 1

[Star Strider]

I generally don’t design filters with such specificity, and so usually use relatively straightforward command-line functions. To design a filter such as yours, I would use the designfilt function (introduced in R2014a). It has a number of arguments, as well as ‘Filter Designer’, a GUI interface that can make this easier. (If you use Filter Designer, click on the ‘File’ tab and then click on ‘Generate MATLAB Code’ to use the filter code wtih the filtfilt function and your signal.) Use filtfilt to do the actual filtering regardless of how you design and code the filter.

[David W Purcell]

Thanks for your help! I am making good progress using designfilt.

Would you agree with the following choices for filter over filtfilt and my choice of filter order?

I normally use filtfilt, but in this case, I think I should use filter because I want the phase effects that the physical GRASS LP511 amp introduces. Their manual says the highpass with 3 Hz cutoff has a fall time constant of 30 ms. Four time constants should get to within 1.8% of a steady state value. My sample rate is 8000 Samples/sec. FIR group delay is 1/2 the number of filter coefficients, so I thought I would use a filter order of:

0.030 sec per Tau * 4 Tau * 8000 Samp/sec * 2 (convert desired filter group delay to filter order) = 1920 filter order or number of filter coefficients

Please let me know if I am in error.

Thank you, DP

[Star Strider]

My pleasure!

The phase delay depends on the type of filter design it uses. (An IIR Bessel design would have no phase distortion in the passband, so in that instance, filtfilt would be preferable. Any other design would have phase distortion.)

I see no problems with your approach. A filter order of 1920 is a long filter, however it may be necessary to get the result you want. I would use the shortest filter possible that would give the desired characteristics. That will obviously require experimentation.

Since you’re designing with group delay as a consideration, also consider the grpdelay function to see what help it may be able to provide in your design process. It returns the values at various frequencies for the supplied filter.

[David W Purcell]

Thank you again. You have helped me realise that I wasn't thinking properly about what my desired phase response charactericsts are. I would like to copy the analogue GRASS LP511 high pass filter characteristics as closely as possible. With your guidance, I am able to achieve the magnitude characteristic quite closely:

FIR_GRASS3Hz=designfilt('arbmagfir', 'FilterOrder', filterOrder, 'Frequencies', manualFreq, 'Amplitudes', manualNormalizedAmplitudes, 'SampleRate', 8000, 'DesignMethod', 'freqsamp');

The GRASS phase response however, is probably like two passive RC filters in series. Each RC section would be expected to have zero phase shift in the pass band, +45deg at the cutoff frequency, and +90deg approaching zero Hertz:

https://www.electronics-tutorials.ws/filter/filter_3.html

I believe that I cannot achieve this phase response using an FIR filter in MATLAB since the constant delay of an FIR filter (half the filter order) causes a linear phase shift for frequency f Hertz at sample rate 8000 Samples/sec of:

phase shift (radians) = f (cycles/s) * -filterOrder/2 (Samples) / 8000 (Samples/sec ) * 2pi (rad/cycle)

Is it possible to simultaneously match the GRASS filter's magnitude and phase respones using an IIR filter?

It may be that the FIR using filtfilt is the closest that I can get to the real GRASS EEG amplifier's high pass. I'd have to accept that the phase response is not being matched.

I'd appreciate any feedback you might be able to provide! Thanks so much for your help, DP

Aside: Other approaches

I believe I could highpass filter using a Fourier transform, applying my gains and phase shifts in the frequency domain, and then inverse Fourier transform. However, my EEG signals are relatively broad band and vary significantly over time. The effects of finite window duration for the Fourier transform would be worse than the FIR choice above.

I also don't want to go to the complexity of a full RC system simulation unless its simpler than it looks!:

https://www.mathworks.com/help/simscape/ug/rc-circuit-in-simulink-and-simscape.html

[Star Strider]

My pleasure!

If you have the Bode plot of the filter, you can use the System Identification Toolbox functions (beginning with the idfrd function), or the Signial Procesing Toolbox invfreqs function, since this appears to be an analogue filter, to estimate its characteristics. (A step or impulse response would do as well, requiring a slightly different approach, however you may not have that information.) That can give you the essential characteristics with respect to pole and zero locations from the transfer function. After that, creating the filter from that information in terms of component values is possible, however definitely not trivial.

These techniques can help you decide if a FIR or IIR filter is appropriate. An IIR design producing a linear phase characteristic in the passband would be a Bessel filter (the besself function).

[David W Purcell]

Again, your help has pointed me in the direction of achieving what I needed! Thank you! I will accept your answer as it answers the FIR question that I originally posted.

With your guidance, I learned what my actual need was and how to solve it. My analogue GRASS filter's measured amplitude response very closely matches an analogue second order highpass system. The function invfreqs that you mentioned above would allow me to use meaurements of amplitude and phase to estimate the s-domain filter coefficients and then I could convert those to z-domain.

mag=labAmp; %measured system amplitude response; using variable name from MATLAB documentation

phase=labPhaseDeg*pi/180; %measured system phase response; deg to rad

w=labFreq*2*pi; %Hz to rad/s; measurement frequencies

h=mag.*exp(1j*phase); %complex response (the transfer function)

[bdat,adat] = invfreqs(h,w,2,2,[],4); %s-domain filter coefficients

However, since the model of a second order high pass system was so good, I decided to directly use that to compute analogue s-domain filter coefficients and then convert those to digital z-domain coefficients for use with the filter function (to obtain the phase response too). I include some code bits below in case they help future readers.

b=[ 1 0 0 ]; w0=3*2*pi; %3Hz cutoff a=[ 1 2*w0 w0^2 ]; %highpass filter coef. with 3 Hz cutoff

[h,w]=freqs(b,a, measuredFreq*2*pi); %model frequency response

[bd,ad]=bilinear(b,a,fs); %convert analogue filter coef. to digital

Numerator coef. in bd, and denominator coef. in ad, can be used with filter(bd, ad, y) to reproduce the filtering characteristics of my real world GRASS LP511 Amp's highpass filter with cutoff 3 Hz.

[Star Strider]

As always, my pleasure!

Thank you for the follow-up.