CODEWITHSUNDEEP
matlabsignal-processingfftaccelerometer
sebap123
sebap123

Reputation: 2685

Accelerometer with FFT - strange output

After reading a lot of research and works on subject I still have got problems applying FFT to my accelerometer data. Most of my code is taken from official MATLAB example: FFT for one dimension. After more reading I've found this question: FFT and accelerometer data: why am I getting this output? where there was suggestion to use windowing. So after some more reading I've added hamming window to my code.

My data looks like that on plot: enter image description here

And this is the code that I am using for FFT:

fs = 1/0.02; %0.02 comes from picking sample each 20ms
m = size(data,1);
w = hanning(m);
yw = w.*data;
n = pow2(nextpow2(yw));
y = fft(yw,size(n,1));
f = (0:size(n,1)-1)*(fs/size(n,1));
power = y.*conj(y)/size(n,1);
figure
plot(f,power)

The problem is that my plot from this code looks like that: enter image description here

Can someone tell me what is wrong with my code? To be honest I'd excepted it would look better (something like this:https://i.sstatic.net/gNdag.jpg) so that's why I am asking this question.

EDIT: My data can be found here: https://dl.dropboxusercontent.com/u/58774274/exp.txt

Upvotes: 1

Views: 4073

Answers (3)

Rash
Rash

Reputation: 4336

Your Fs is 50 so the highest frequency in your data can be Fs/2 = 25Hz.

See if this code helps.

fid = fopen('1.txt','r');
C = textscan(fid, '%f');
fclose(fid);
data = C{1};
fs = 50;
m = length(data);
nfft = 2^nextpow2(m);
y = fft(data,nfft)/m;
f = fs/2 * linspace(0,1,nfft/2+1);
power = abs(y);
subplot(211)
plot(f,power(1:nfft/2+1)) 
t = (0 : m-1)/fs;
s0 = .8*fs : 3.2*fs;  % .8 sec to 3.2 sec
p(s0) = .5*cos(2*pi*3.3*t(s0)+.25*pi); 
p = p + mean(data); 
subplot(212)
plot(t,data);hold on
plot(t,p,'r')

enter image description here

This is your data in frequency domain.

There is a peak at 3.3 Hz.

As a proof I plotted a sinusoidal with frequency of 3.3 Hz along with your data,

As you can see, it totally matches your data.

Upvotes: 2

hotpaw2
hotpaw2

Reputation: 70693

Your plots will look better if you first remove the DC offset (subtract the average of all samples from each point before computing the FFT), and then plot only half the FFT result data points (N/2) or less (since the upper half of an FFT result is just a conjugate mirror of the first half for real data input).

Upvotes: 2

lmNt
lmNt

Reputation: 3892

What your reference plot shows is just the magnitude on a logarithmic scale (in dB). Also it seems that the DC offset is removed. So you end up with something like this if you adapt it in your code.

data = data-mean(data); % DC removal
fs = 50;
m = length(data);
nfft = 2^nextpow2(m);
y = fft(data,nfft)/m;
f = fs/2 * linspace(0,1,nfft/2+1);
power = abs(y);
plot(f,10*log10(power(1:nfft/2+1))); % plot log magnitude
ylim([-30 0]) %limit axis

Looks at least similar to me.

enter image description here

If you additionally take the absolute square you are basically estimating the power spectral density of your signal. So meaning you have to change power = abs(y); to power = abs(y.^2); and then obviously have to adapt your ylim([-30 0]) to something lower like -50.

Upvotes: 1

Related Questions