butterworth filter cutoff frequency calculation

Hi All, I am not a electrical engineer so my question may sound stupid. Sorry about that. I have a some noisy data and I want to filter my data using a butterworth filter of order 2. I am having trouble finding the cutoff frequency for the filter design. How do I calculate cutoff frequency? I am sampling my data every 10 seconds. Thnanks.

Answer 1

[Star Strider]

First, a Butterworth (or any) filter of order 2 is not likely to work as well as you want it to. You may find a Chebyshev Type II design more appropriate, especially with your extremely low sampling frequency.

Second, the best way to design a band-selective filter is to take the fft (link) of your signal to determine where the signal frequencies are and where the noise frequencies are. Use that information to design your filter. Note that if you have broadband noise, a frequency-selective filter will not eliminate all of the noise. You might want to eliminate baseline offset and low-frequency baseline variation as well, the reason I usually use a bandpass design.

The designfilt function can make filter design straightforward. Another acceptable way to design a Butterworth filter is this bandpass prototype:

Fs = 0.1;                                                   % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency
Wp = [0.002  0.010]/Fn;                                     % Specify Bandpass Filter
Ws = [0.001  0.015]/Fn;                                     % Stopband (normalised)
Rp = 10;                                                    % Passband Ripple (dB)
Rs = 30;                                                    % Stopband Ripple (dB)
[n,Wn] = buttord(Wp, Ws, Rp, Rs);
[z,p,k] = butter(n,Wn);                                     % ZPK 
[SOS,G] = zp2sos(z,p,k);                                    % Convert To SOS for Stability

figure(2)
freqz(SOS, 2^16, Fs)

To design a Chebyshev Type II filter, replace buttord with cheb2ord and butter with cheby2. The arguments are slightly different, so keep that in mind. The rest of the code does not change.

[masih]

Just a follow up question this please: what is the criteria to find Wp Ws Rp and Rs? I have attached a sample of my data. Thanks

[masih]

here it is

[Star Strider]

The values for ‘Ws’ and ‘Wp’ derive from your data. You have to decide those, based on the Fourier transform of your signal. The values for ‘Rp’ and ‘Rs’ are also empirical, so they have to design a stable filter that does what you want. (They are essentially irrelevant in a Butterworth design, but are relevant in Cheybshev and other designs.) I do not know what you want from your filter, so I leave the exact values of the passband and stopband frequencies for you to experiment with.

The Code —

[d,s] = xlsread('masih data_excel1.xls');
t = datenum(s(2:end,1), 'yyyy-mm-dd HH:MM:SS');

figure(1)
plot(t, d)
grid
datetick('x', 'HH:MM:SS', 'keepticks')

Ts = mean(diff(t))*(24*60*60);                          % Sampling Time (sec)
Fs = 1/Ts;                                              % Sampling Frequency (Hz)
Fn = Fs/2;                                              % Nyquist Frequency (Hz)
L = length(t);
dsm = d-mean(d);                                        % Subtract Mean
FTd = fft(dsm)/L;                                       % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                     % Frequency Vector
Iv = 1:length(Fv);                                      % Index Vector

figure(2)
plot(Fv, 2*abs(FTd(Iv,:)))
grid
axis([0  0.01    ylim])

Wp = [0.0002  0.0070]/Fn;                               % Specify Bandpass Filter
Ws = [0.0001  0.0075]/Fn;                               % Stopband (normalised)
Rp = 10;                                                % Passband Ripple (dB)
Rs = 30;                                                % Stopband Ripple (dB)
[n,Wn] = buttord(Wp, Ws, Rp, Rs);
[z,p,k] = butter(n,Wn);                                 % ZPK 
[SOS,G] = zp2sos(z,p,k);                                % Convert To SOS for Stability

figure(3)
freqz(SOS, 2^16, Fs)                                    % Filter Bode Plot

d_filt = filtfilt(SOS, G, d);                           % Filter Signal

figure(4)
plot(t, d_filt)
grid
datetick('x', 'HH:MM:SS', 'keepticks')

[masih]

Thanks so much. This is really helpful. Can you please tell me how you came up with [0.0002 0.0070] and [0.0001 0.0075] for Wp and Ws? Thanks a lot.

[Star Strider]

As always, my pleasure.

I came up with the values for ‘Wp’ by looking at the fft plot. I chose the ‘Ws’ values to be just ‘outside’ (less than and greater than, respectively) the ‘Wp’ values.

I am not certain what you want from your data. Experiment with the lower (first) values for ‘Wp’ and ‘Ws’ to get what you want. I set the lower values to eliminate the d-c (constant) offset and low-frequency baseline variation. (The higher values filter out some of the high-frequency noise.)

If you want to keep the baseline offset and low-frequency variation, it may be best to design a lowpass filter that just eliminates the high-frequency noise, instead of the bandpass filter I designed. See the documentation for the butter function for information on how to do that. The only change you would need to make is for ‘Wp’ and ‘Ws’ to be scalars (rather than vectors), perhaps using the higher values I used in my bandpass design to start with. The default filter design is lowpass, so you would not need to specify that.

Experiment with my code to get the result you want. I will help as I can.

[masih]

Thanks a lot. I guess my question is mainly how we can find the value of Wp by looking at the fft plot. Can you please show me on the fft plot how you got that for my data? Appreciate it.

[Star Strider]

My pleasure.

Your data have some large, low-frequency peaks and higher-frequency noise. I actually picked ‘Wp’ arbitrarily, based on what I believed were the ‘significant’ frequencies. Signal processing of data is essentially just ‘cut-and-try’ until you get the result you want.

Look at the Fourier transformed data (in my code, figure(2)), and see how a cutoff frequency of 0.0005 would work with your data in a lowpass filter design:

Wp = 0.00050/Fn;
Ws = 0.00060/Fn;

(Nothing else in my code changes.) Then adjust the passband and stopband frequencies up and down to see how the result changes. When you are happy with the result, you have defined your frequencies-of-interest in your data.

Using findpeaks:

[pks, frqs] = findpeaks(2*abs(FTd(Iv,1)), Fv, 'NPeaks',3)

your data have three ‘significant’ peaks (in my opinion), at frequencies of 0.0001001, 0.0002002, 0.0004004. I used those to set the passband frequency, and slightly increased it to define the stopband frequency. What you choose depends on what you want from your data.

[Ahmed Hassan]

Hey,

Could you please help me !!

i was trying to understand your code but I did not get how to enter my values.

for example If I want to know the cutoff frequency for my analoge Butterworth filter for the following given data:

max frequency i will get in the passband is 9 kHz and i want to attenuate frequency higher than this value.

sampling frequency 20 kHz of even 40 kHz

stopband Ripple not less than 55 db

what could be my code to know the cutoff frequency ?

[Star Strider]

@Ahmed Hassan — Post this as a new Question, and include more information.

I will delete this Comment in a few hours.

[Daniel M]

Hi Star Strider, I'm trying to follow along, but when I run your code, I get the following error:

Warning: Usage of DATENUM with empty date character vectors or empty strings is
not supported.
Results may change in future versions. 
> In datenum (line 95)
  In untitled48 (line 2) 
Error using plot
Vectors must be the same length.

Error in untitled (line 4)
plot(t, d)

t is empty because datenum() failed.

It seems to be an issue with xlsread on Linux, because

whos d s
  Name         Size            Bytes  Class     Attributes

  d         1999x3             47976  double              
  s            1x3               396  cell 

On Windows it works properly. So nevermind then, I'll just do it there.