Prev: error in simulink
Next: how to plot gamma distribution
From: Matthew on 23 Oct 2009 10:49
I want to create a lowpass butterworth filter for some data that I have and from what I've read online I can't figure out how to actually implement the filter. I have time and acceleration data (so two vectors, say, Time and Accel) collected at 500kHz and I'd like to do a lowpass filter with a cutoff frequency of 20kHz with no phase shift. Not sure about the order yet, but will assume 5th order right now.

Everything I found suggests that I should use [b, a] = butter(5,500000,'low'). But I'm not sure what to do with [b, a], whether I should apply it to the time and acceleration vectors, what the output would be, and how I could plot this signal on top of my original vectors.

Any help would be greatly appreciated.
From: Wayne King on 23 Oct 2009 11:08
"Matthew " <mdford86(a)gmail.com> wrote in message <hbsfpg$c3g$1(a)fred.mathworks.com>...
> I want to create a lowpass butterworth filter for some data that I have and from what I've read online I can't figure out how to actually implement the filter. I have time and acceleration data (so two vectors, say, Time and Accel) collected at 500kHz and I'd like to do a lowpass filter with a cutoff frequency of 20kHz with no phase shift. Not sure about the order yet, but will assume 5th order right now.
>
> Everything I found suggests that I should use [b, a] = butter(5,500000,'low'). But I'm not sure what to do with [b, a], whether I should apply it to the time and acceleration vectors, what the output would be, and how I could plot this signal on top of my original vectors.
>
> Any help would be greatly appreciated.
Hi Matthew, the outputs b and a are the numerator and denominator coefficients of the IIR filter you design with butter(). If you have collected acceleration data as a function of time, then you want to apply the filter to the acceleration data. It is not possible to design a realizable filter with zero phase. However, you can do "zero phase" filtering in Matlab with the filtfilt() function. See

>>doc filtfilt
% you should read the documentation

The relevant syntax is
output = filtfilt(b,a,input_data); % the input_data are your acceleration %measurments

However, you should be aware that butter() works on normalized frequency. Your syntax for butter() '[b, a] = butter(5,500000,'low')' can't be correct. The 2nd input should be your cutoff frequency (not the sampling frequency) and the cutoff frequency has to be in the interval [0,1] with 1 corresponding to pi radians/sample. You have to convert your cutoff frequency of 20 kHz to normalized frequency:

(2*pi*2e4)/5e5 %assuming your sampling frequency is 500 kHz as you say.

Hope that helps,
wayne
From: Matthew on 23 Oct 2009 12:05
Thank you for the help. One last question: If I create a filter using the fdatool and save it as an M-file, how do I apply this filter to a vector in my workspace? For example, if I create a butterworth filter and call it lowpassbutter, I get this M-file:

function Hd = lowpassbutter
%LOWPASSBUTTER Returns a discrete-time filter object.

%
% M-File generated by MATLAB(R) 7.8 and the Signal Processing Toolbox 6.11.
%
% Generated on: 22-Oct-2009 17:26:32
%

% Butterworth Lowpass filter designed using the BUTTER function.

% All frequency values are in Hz.
Fs = 500000; % Sampling Frequency

N = 15; % Order
Fc = 20000; % Cutoff Frequency

% Calculate the zpk values using the BUTTER function.
[z, p, k] = butter(N, Fc/(Fs/2));

% To avoid round-off errors, do not use the transfer function. Instead
% get the zpk representation and convert it to second-order sections.
[sos_var,g] = zp2sos(z, p, k);
Hd = dfilt.df2sos(sos_var, g);

% [EOF]


At this point, how do I call the filter? I tried outputvector = filter(lowpassbutter, inputvector) but I get vector full of NaN. Any help would be greatly appreciated.
From: Wayne King on 23 Oct 2009 12:26
"Matthew " <mdford86(a)gmail.com> wrote in message <hbsk80$12a$1(a)fred.mathworks.com>...
> Thank you for the help. One last question: If I create a filter using the fdatool and save it as an M-file, how do I apply this filter to a vector in my workspace? For example, if I create a butterworth filter and call it lowpassbutter, I get this M-file:
>
> function Hd = lowpassbutter
> %LOWPASSBUTTER Returns a discrete-time filter object.
>
> %
> % M-File generated by MATLAB(R) 7.8 and the Signal Processing Toolbox 6.11.
> %
> % Generated on: 22-Oct-2009 17:26:32
> %
>
> % Butterworth Lowpass filter designed using the BUTTER function.
>
> % All frequency values are in Hz.
> Fs = 500000; % Sampling Frequency
>
> N = 15; % Order
> Fc = 20000; % Cutoff Frequency
>
> % Calculate the zpk values using the BUTTER function.
> [z, p, k] = butter(N, Fc/(Fs/2));
>
> % To avoid round-off errors, do not use the transfer function. Instead
> % get the zpk representation and convert it to second-order sections.
> [sos_var,g] = zp2sos(z, p, k);
> Hd = dfilt.df2sos(sos_var, g);
>
> % [EOF]
>
>
> At this point, how do I call the filter? I tried outputvector = filter(lowpassbutter, inputvector) but I get vector full of NaN. Any help would be greatly appreciated.

Your function lowpassbutter.m should return a dfilt filter object called Hd. Make sure that lowpassbutter.m is in your Matlab path and run
>>Hd = lowpassbutter;
Then you should be able to type
>>whos Hd
and see that there is a dfilt object, specifically a dfilt.df2sos object in your Matlab workspace. You can use this dfilt object as an input to filter()

output = filter(Hd,input_data); %filter your input data with the filter object

However, you should realize that filter() does NOT implement zero-phase filtering as you indicated you wanted in your original post. If you're okay with this, fine. If you need to use filtfilt(), then you will need to have the filter coefficients because filtfilt() does not currently accept filter objects as inputs.

wayne
 | 
Pages: 1
Prev: error in simulink
Next: how to plot gamma distribution