2
votes

I'm having difficulty in calculating the Butterworth coefficients with the 'buttord' and 'butter' functions. My objective is to filter out noise from a time series I constructed. The time series has a red noise component and a sinuosoid signal with frequency 0.3 Hz. The sampling frequency of the time series is 10 Hz.

Following the documentation on 'buttord' http://www.mathworks.com/help/signal/ref/buttord.html I calculated [n, Wn] for the specifications (followed example 1 of the documentation):

Wp = 0.33/(sfreq/2); Ws = 0.37/(sfreq/2);
passripp = 0.1; stopatten = 40;
[n,Wn] = buttord(Wp,Ws,passripp,stopatten);
[b,a] = butter(n,Wn);
y_butter = filter(b,a,timeseries(:,2));

Plotting y_butter with time just gives me zero everywhere!

I attempted to use 'freqz' to examine the frequency response of the filter (using 512 samples):

freqz(b,a,512,sfreq)

whose plot shows that the transition band is between 1 and 4 Hz!

My understanding behind the filter is:

  • signal at 0.3 Hz
  • noise at >> 0.3 Hz
  • pass everything from 0 to 0.33 Hz
  • attenuate everything from 0.36 Hz beyond

Your help will be much appreciated!

Data can be downloaded here: http://dl.dropbox.com/u/1918592/detrendedTS.mat Column 1 of 'ts' is time variable, column 2 is data variable

I detrended a linear fit (Matlab 'detrend') to remove some of the walk-away red noise behavior.

1

1 Answers

3
votes

The code that you posted generates a 57th (!!!!) order filter to perform the operation that you want to do. The freqz command does indeed show that this filter is valid, but attempting to implement such a high order filter directly without first converting it to 2nd order sections using tf2soswill undoubtedly result in serious numerical errors. This is probably why you just see zeros. If you convert the filter coefficients into 2nd order sections and then cascade the filters, you should actually get the filter output observed using the freqz command.

The reason you get such a high filter order in the first place is because you have placed the start and stop bands very closely. The use of the buttord function is completely unnecessary for a simple low-pass Butterworth filter. Just use ...

[b,a] = butter(n,cut/(sfreq/2));

... just pick n to give you the rolloff you want (6,12,18 dB etc . . ). You'll only run into problems once n becomes about 10 when using double precision data representation, but you'll probably find a relatively small filter order will do the job you want.