Signal Processing Help

I am a post-bac research assistant. My background is in physiology. I am new to signal processing. I have been asked to write a script to reduce the noise in a data set and run it on several files.

The process as I understand it... 1. Convert data from time to frequency domain. 2. Filter data with Butterworth filter. 3. Convert data from frequency to time domain.

Example of my code...

N=length(Data);
X = abs(fft(Data,N));
F = [0:N-1];
plot(F,X)
xlabel('Frequency / fs')
ylabel('Amplitude')

[B,A] = butter(2, 0.5);

figure
Xf=filter(B,A,X);
plot(F,Xf,'b-');
title('Filterd Frequency Data')
xlabel('Freq (Hz)')
ylabel('Amplitude')

figure
Nf=length(Xf);
Y = abs(ifft(Xf,Nf));
F = [0:Nf-1];
plot(F,X);

I have tried watching and reading tutorials online but am very lost. Would love help understanding how to use fft and butter. Any tutorials, guides, example code or suggestions would be great! I have the signal processing toolbox. Thank you for your help!

Some questions that come to mind...

 - Am I using fft correctly?
 - Do I need to use abs or fftshift with fft?
 - How do I understand how to set the cutoff frequency?
 - How do I determine the order of the Butterworth filter?
 - Am I using filter correctly?
 - Am I using ifft correctly?
 - Do I need to use abs or fftshift with ifft?

Answer 1

[Wayne King]

Generally the stricter your filter requirements the higher the order.

The relationship between any linear time invariant operator, like your filter, and the Fourier transform is this one: filtering an input signal with your filter is equivalent to multiplying the Fourier transform of your signal with the Fourier transform of the filter.

Knowing what filter specifications you need is up to you, you have to know what the spectrum of your signal is, what frequency range are you interested in, what frequencies would you like to reject. How much noise is there in the band of frequencies you would like to reject. Those answers are specific to your application.

Answer 2

[Wayne King]

Hi, You do not want to use filter() on the frequency domain representation of the data. Use filter() directly on the time domain signal.

output = filter(B,A,data);

It's hard to say whether you have designed your filter correctly because you do not give the relevant details such as the sampling frequency and what frequency you want to filter below (cutoff) frequency.

Your filter:

[B,A] = butter(2, 0.5);

is a very relaxed one, it only attenuates the input by 3 dB at a frequency 1/4 of the sampling rate. Let's assume your data is sampled at 1000 Hz, to see what your filter's magnitude response looks like:

fvtool(B,A,'Fs',1000);

It's basically passing everything from [0,250] Hz unattenuated. What specifications do you want?

[Daniel]

Thanks! So I do not need to use ifft? I can filter both ‘data’ and the associated ‘time’ variables...

output1 = filter(B,A,data);

output2 = filter(B,A,time);

Ideally, my data is recorded 200 times per second (200 Hz?). However, the actual number is usually lower depending on the number of electrodes being recorded. I currently have zeros as place holders for the missing values. Do I need the place holders or should I delete them?

Is my first figure correctly plotting frequency? I saw it done several different ways online and find it confusing. If it is correct, a stem plot of my first figure shows a few large protruding vertical lines on the left hand side, with small consistent vertical lines comprising the rest of the plot. I am assuming that I want to eliminate the larger vertical lines on the left hand side. Can I simply use the X (frequency) value from the stem plot as my cut off value?

[b,a]=butter(n,Wn)

What does ‘n’ mean? How do I determine ‘n’?

Would ‘Wn’ be 'Cut off frequency'/500?

Thank you for your help!

[Daniel]

Updated code…

%data = x-velocity;

%time = time of velocity recording (200 recordings per second);

%state= was the value recorded at give time (1=Yes & 0=No);

% Converted data from time to frequency domain

N=length(data);

X = abs(fft(data,N));

F = [0:N-1]/N;

plot(F,X)

xlabel('Frequency / fs')

ylabel('Amplitude')

%Apply Butterworth filter… tried/guessed many different values…

[B,A] = butter(7,.08);

% Filtered values

data_out = filter(B,A,data);

time_out = filter(B,A,time);

state_out = filter(B,A,state);

% Delete values if not recorded (state==0)

data_out=data_out(state_out>0.5);

time_out=time_out(state_out>0.5);

% Plot data

plot(time_out,data_out);

title('Filterd Velocity Data')

xlabel('Time')

ylabel('Velocity')

Still not there but hopefully getting closer...

Answer 3

[Wayne King]

Hi, you do not need to use ifft(). You do not need to filter the time vector, just the data.

output = filter(B,A,data);

butter() and many other MATLAB filter design functions always use normalized frequency from [0 1] where every frequency in that interval is really 1/10*pi, 1/3*pi, 1/2*pi, 1*pi. In other words whatever number you put in for your cutoff frequency is really that number times pi.

This frequency is in radians/sample. Usually people want to specify their filter design in Hz. You can convert Hz to normalized frequency if you know the sampling frequency.

Wn = (2*f)/Fs;

where Fs is the sampling frequency.

Ex:

Assume the data is sampled at 1 kHz and you want the cutoff frequency at 200 Hz.

Wn = (2*200)/1000;
[B,A] = butter(10,Wn,'low');
fvtool(B,A,'Fs',1000);

Of course you have to determine what order and cutoff frequency make sense for you.

You can use fdesign.lowpass to design your filter directly in Hz.

d = fdesign.lowpass('N,F3db',10,200,1000);
Hd = design(d,'butter');
output = filter(Hd,data);

[Daniel]

Now I'm starting to break some ground!

I am able to filter my data (in the time domain) with a butter filter and smooth the line. I have changed the order and frequency of the butter filter (guess and check) and seen improvement.

Another lab use to process the kinematics of our experiments. We are processing our own kinematics for the last experiment. The other lab suggested we use a fourier transform and a butterworth filter to reduce signal noise.

What is the connection between fft and butter?

I can achieve okay results with 6th order and f=5... but I would love a more scientific way of deciding filter criteria than guess and check!

How do I determine the order and frequency cutoff of my butter filter?

Thank you for your help!