Help filtering a signal using fft

64 次查看（过去 30 天）
simon
simon 2016-10-9
评论： Star Strider 2016-10-10
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!
  2 个评论
John BG
John BG 2016-10-9
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
simon 2016-10-9
Good point, should have thought of it first... Thanks for the heads up, and it is attached now!

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2016-10-9
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
  2 个评论
simon
simon 2016-10-10
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
Star Strider 2016-10-10
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;

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Estimation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by