Help filtering a signal using fft

Hey,

I'm trying to filter a signal with the use of the fft. I would like to keep the lower frequencies and cut the higher frequencies.

I have given the code below.

I assume that I'm wrong in defining my filter variable... Right now, my FFT has length ~8000, and I wanted to keep the first 300 samples of that fft, and cut off the rest.

clear all;
close all;
load('dataTP.mat');

L=length(dataTP.data);

%create data vector, the one I want to filter
for i=1:length(dataTP.data)
    data(i)=sqrt((dataTP.data(i,1))^2+(dataTP.data(i,2))^2+(dataTP.data(i,3))^2);
end

data = data-mean(data);

%create filter for use in fftfilt
filter=zeros(1,length(data));
filter(1:300)=1;

filtered_data=fftfilt(filter, data);

%create fft to visualize in graph
data_fft=fft(data);
P2 = abs(data_fft/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);

f = 200*(0:(L/2))/L;

%plot fft visualization
figure(1)
plot(f,P1)

%plot raw data
figure(2)
plot(data);

%plot filtered data vs raw data
figure(3)
plot(F)
hold on
plot(data, 'r')
hold off

My first attempt was made like this:

data_fft=fft(data);
filter=zeros(1,length(data));
filter(1:300)=1;  

filtered_fft=data_fft.*filter;
filtered_data=ifft(filtered_fft);

but that didn't work too well either... Hope you can help!

I have attached the data File.

Thanks!

[John BG]

could please make available the dataTP1.mat file, or a portion of it so the readers of your question can reproduce your set-up?

[simon]

Good point, should have thought of it first... Thanks for the heads up, and it is attached now!

Answer 1

[Star Strider]

Ideally, you should also have a time vector. This allows you to define specific frequencies rather than referring to elements of the Fourier transform vector. I created a time vector, so you can substitute it for yours or leave my code as it is.

See if this does what you want:

d = load('simon DataTP.mat');
DataRaw = d.dataTP.data;
data = sqrt(sum(DataRaw.^2,2));                                     % Euclidean

Ts = 1;                                                             % Sampling Interval
Fs = 1/Ts;                                                          % Sampling Frequency
Fn = Fs/2;                                                          % Nyquist Frequency
L = size(data,1);                                                   % Length Of ‘data’ Vector
t = 1:L*Ts;                                                         % Time Vector

FTdata = fft(data)/L;                                               % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                                 % Frequency Vector (One-Sided FFT)
Iv = 1:length (Fv);                                                 % Index Vector
Iv300 = 1:300;                                                      % ‘First 300’ FFT Frequencies
F300 = Fv(300);                                                     % Frequency Corresponding To 300th Element

figure(1)
plot(Fv, abs(FTdata(Iv))*2)
grid

figure(2)
plot(Fv(Iv300), abs(FTdata(Iv300))*2)
grid

FLen = 48;                                                          % Discrete Filter Order
b_filt = fir1(FLen, F300/Fn);                                       % Design FIR Filter

figure(3)
freqz(b_filt, 1, 4096, Fs)

data_filtered = fftfilt(b_filt, data);

figure(4)
subplot(2,1,1)
plot(t, data)
title('Original Data')
grid
subplot(2,1,2)
plot(t, data_filtered)
title('Filtered Data')
grid

[simon]

This works great, thanks!

However, our sampling frequency (on the sensors) is done at 200 Hz, and when I enter that, it screws up the time vector... iI did that just by replacing the Ts=1 with Ts=1/200. What am I doing wrong?

Thanks again though!

[Star Strider]

My pleasure!

Overnight, my primary HP Win 8.1 MATLAB computer died (it needs yet another power management chip, as it seems to need every couple years), and this Win 10 computer (running R2015b) won't let me load .mat files with either Firefox or IE. (My opinion of Win 10 is not fit for an open forum.) So I can’t load your file to test it. I'll access Answers on another machine or on my phone to see if I can download it and then transfer it via Bluetooth later. I’m tired of having to do work-arounds for Win 10’s myriad deficiencies.

See if this works for ‘t’:

t = linspace(0, 1, L)*Ts;