Processing Accelerometer Vibration Data With FFT

Question

For a project I am attempting to use an accelerometer to measure vibration in an RC aircraft and determine frequency from the result. To do this I am using an MPU-6000 accelerometer sampling @ 1000Hz. Before I actually deploy this system on an aircraft, I have taped the sensor to a speaker and used a frequency generator to generate a 100Hz tone for 30 seconds. This should be well within the accelerometer's range:

1000Hz / 2 = 500Hz maximum by Nyquist rule, much higher than 100Hz

My amplitude data in the time domain is bellow (generated w/ GNU Octave):

Zoomed in:

That looks like a ton of data to me but nothing seems to be obviously in error. Bellow is the data after performing a simple periodogram in GNU Octave with the command periodogram(DATA,[],[],1000):

A peak can be seen at 100Hz but there are all of those other spikes that I can't explain. I have tried to detrend my data using the GNU Octave command detrend(DATA,'linear') but that had no effect on removing the stray peaks.

I get the feeling that I should be doing some filtering but I am new to this and really don't know what I am doing. How can I remove all of the stray peaks in the periodogram? Thanks in advance!

Data file on Google Drive: http://goo.gl/w3Brol

EDIT:

Here is the periodogram for the unit sitting on the speaker with nothing playing for 30 seconds:

There doesn't appear to be any interference. The speaker was on, just not playing anything. This is making me think that my speaker is generating harmonics when I supply a sine wave. Could that be? As always, thanks for the help!

Silence data file on Google Drive: http://goo.gl/bNbnDA

EDIT #2:

Success! After isolating my speaker from the table, taping the device to the back of the speaker (I bet glue would help further), and lowering the volume I was able to get the result below! You guys were spot on, thanks for all your help!

Answer 1

Back when I used to design accelerometers, our biggest test issues were mechanical mounting stiffness, and the aforementioned ground loops.

1. Eliminate ground issues - The 50Hz/60Hz ground loops can be eliminated by taking the data with no power line attached. Unplug the laptops, so it runs off batteries. If you must have line-powered instrument(sig gen, DAQ, amplifier), power it from a UPS(unplugged). You can rule out ground loops if ALL the harmonics move when you sweep the frequency down. If there is a ground loop, it will stay put at 50/60Hz and the harmonics of that.
   1a.Try inputting odd signals of 53Hz or 59Hz, something whose harmonics don't coincide with 50Hz 60Hz, or your sampling freq. If you might be aliasing, change your sample freq, and see if the harmonics move.
2. Mounting. I didn't download the data, but the amplitude dropped during the time domain plot. If the accel loosened up, and is rocking, even microscopically, it can create harmonics as it impacts on the left/right side, every up/down cycle. These can mix back down to subharmonics, like 50Hz. I'm guessing the accel did not stay solidly attached to the speaker cone. As a mount, tape usually won't cut it at any frequency. Here's a link illustrating how even metal studs aren't enough to accurately reproduce high frequencies. http://www.mmf.de/accelerometer_mounting.htm.
   2a. There is also mass loading to consider, if you are adding the accel to the edge of a flimsy paper speaker cone, which can barely drive its own weight. At 100Hz and above, you may get a more solid mount on the back of the speaker magnet. Adhering with thin layer of beeswax, or superglue, might get you to 1KHz, without requiring a machined mount. While back-mounting will be admittedly less signal in g, it will be a cleaner one, from a distortion perspective(assuming a distortion free source).

If its not convenient to retake the data, try a rolling STFT, or spectrogram. If the harmonics jump up part way through duration of the measurement, I'd look hard at the mount.

Answer 2

After taking a closer look at your results, here is what I think.

Because you are playing back sinusoid through your speaker you are also observing it's properties. It has some significant harmonic distortions appearing as the harmonics of fundamental frequency: $200$, $300$, $400 \mathrm{Hz}$ in your case. That is either because:

• You are driving your speaker at high volume
• It is simply crap for accurate measurements (sorry...)

Obviously you could also get harmonic distortions coming from your sound card, but this is unlikely.

On the other hand what you see at $50 \mathrm{Hz}$ is the mains hum of AC voltage (you are definitely not in USA). This distortion also has it's harmonics at $100$, $150$, $200 \mathrm{Hz}$, etc. So you can see that they do overlap with your test signal.

Nevertheless if you look at higher frequencies, i.e. around $400 \mathrm{Hz}$, you will see two peaks next to each other. Those are two harmonics, one of your test signal and other one of your mains hum. Because there was a slight offset at low frequencies it diverges even more at higher frequencies.

In order to confirm the mains hum you should simply record silence. If it is present, then well - you have a slight problem. If you are using laptop, try to unplug the AC adaptor - it should help. Otherwise you are simply picking up everything with your power supply. You can also try external audio interface. On the contrary, if you are already using some kind of interface with XLR microphone or whatever - try to break the ground loop.

When it comes to speaker distortions try to use other one or maybe use it a slightly lower volume.

Answer 3

This is just a guess, but it could be that you are getting these harmonics because you have inadequate acoustic coupling between the accelerometer and its mount point (i.e. tape isn't secure enough). If you have a way to sense data from a microphone you could rule out the possibility of the speaker being defective.

As far as filtering out the harmonics, I thought this would be a good opportunity to try out a technique I learned recently called Singular Spectrum Analysis. The idea is that you use the singular value decomposition to break apart your signal into a weighted combination of orthogonal basis functions (similar to the fourier basis, though the functions aren't limited to complex exponentials).

The spectral line at 100Hz in your periodogram is clearly the dominant mode. We should be able to extract it with a rank-2 approximation of your signal (rank-2 because a complex exponential at 100Hz has both real and imaginary parts -- so we need 2 real-valued basis functions to represent it).

Below is a matlab script that uses the singular spectrum analysis technique I mentioned above.

clear all
close all

path = '100.csv';
raw = csvread(path);

%% spectral decomposition with svd
J = 100; % number of spectral components
K = length(raw) - J + 1; % length of basis vectors
col = raw(1:J);
row = raw(J:end);
H = hankel(col,row)'; % tall, thin hankel matrix. 
[U,S,V] = svd(H, 'econ');
s = diag(S); % squared singular values of H
stem(s);
grid on
title('singular spectrum');
xlabel('component index');
ylabel('squared singular values');

%% reconstruct signal from the 2 dominant spectral components
Hfilt = U(:,1:2)*S(1:2,1:2)*V(:,1:2)'; % a rank-2 approximation of H 
y = zeros(size(raw)); % this will be the filtered signal
for j=1:J
    sprintf('Processing column %d of %d ...', [j, J])
    for k=1:K
        y(k+j-1) =  y(k+j-1) + 1/J*Hfilt(k,j); % anti-diagonal averaging
    end
end

%% power spectra of original and filtered signals
% I'm using Welch's method, which averages a bunch of lower resolution
% periodograms instead of computing a single high-resolution periodogram.
% This reduces the variance of the power spectrum estimate in exchange
% for lower frequency resolution.
Fs = 1000;
nfft = 1024;
win = blackmanharris(nfft); % this window has really high sidelobe attenuation
noverlap = ceil(.661*nfft); % near-optimal overlap for blackman-harris window
[R, f] = pwelch(raw, win, noverlap, nfft, Fs, 'onesided');
[Ry, f] = pwelch(y, win, noverlap, nfft, Fs, 'onesided');
figure
plot(f, 10*log10(Ry), 'linewidth', 2);
hold on
plot(f, 10*log10(R), 'r');
grid on
legend('filtered', 'original');
xlabel('freq. (Hz)');
ylabel('power (dB)');
title('power spectra');

This plot shows that most of the energy of your signal is concentrated in only 2 modes. These capture your 100Hz complex exponential.

Here's the result of the rank-2 approximation compared to the original signal.