Designing a differentiator filter using cfirpm

Hi,

I would like to use a differentiator filter with response and phase function as shown below as described in this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2235870/

The purpose of the filter is to normalize power spectrum by compensating for the 1/f trend. The author mentioned that their method can be implemented using cfirpm in Matlab but they did not give any further information about the parameters that they used. Any pointers would be appreciated. Thank you.

Maria

Answer 1

[Star Strider]

The moving average filter is easy enough to design:

Fs =  250;                                      % Sampling Frequency (Hz) For EEG
Ts = 1/Fs;                                      % Sampling Interval (sec)
b = [1 -1];                                     % MA Filter Numerator
a = 2*Ts;                                       % MA Filter Denominator

figure
freqz(b, a, 2^14, Fs)                           % Filter Bode Plot

It gives a result similar to the published Bode plot.

A single-band differentiator using cfirpm is:

n = 60;
f = [0.01 0.9];
b = cfirpm(n,f,{@differentiator,Fs})
a = 1;

figure
freqz(b,a,2^14, Fs)EDIT — 

EDIT —

Use either filter or filtfilt to filter your signals with these filters. (They are linear-phase filters, so both will work with them. Use filtfilt if there are problems with the filter initial conditions. or the output from filter is not what you would expect.)

[Maria Wijaya]

Thank you for the answer. I tried using filter and filtfilt on my data but the amplitude becomes very large after filtering. I just changed the sampling rate from your example. Is there anything else that needs to be modified? Code and time series plots below:

Fs = 500;
n = 60;
f = [0.01 0.9];
b = cfirpm(n,f,{@differentiator,Fs});
a = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs =  500;                                    
Ts = 1/Fs;                                      
b = [1 -1];                                     
a = 2*Ts;                                       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eeg = double(EEG.data(2,:,1)); % 1 epoch of EEG data from 1 channel, 1850 data points
eegf = filter(b,a,eeg);
eegfilt = filtfilt(b,a,eeg);

[Star Strider]

My pleasure.

The sampling frequency in the filter designs has to be the sampling frequency of your signals.

I noticed the filter transfer functions amplified the signals, however the paper had similar results.

The way I would deal with that is to change ’b’ for both filters:

a = Ts * 10^3.1;                                % MA Filter Denominator

and for the cfirpm filter:

a = 10^3.1;

These are the magnitude of the maxima of the transfer functions, and are calculated from the maximum amplitude of both being about 62 dB, so . The maximum amplitude of the transfer functions is about 0 dB with these changes, as desired for filters.

That change should stop them from being amplifiers, and produces the correct result when I evaluate them with freqz. Note that they significantly attenuate low frequencies, and (being differentiators) will amplify high-frequency noise. They will also remove baseline offset.

Also, in your code, you are filtering with the moving average filter coefficients. The two filters have similar characteristics, so that may not make a significant difference.

I’ve not used these types of filters with EEG signals (or any other signals, for that matter), generally using elliptic filters to remove band-limited noise and to remove baseline offset and low-frequency baseline wander. It’s certainly an interesting approach to dealing with the noise.

[Star Strider]

If my Answer helped you solve your problem, please Accept it!