Fast Fourier Transform for varying period (Order analysis)

I am trying to do a fast Fourier transformation on accelerometer data from a shaft rotating at varying speed in Python.

What I have done so far:

1: The original plot was in the time domain, and I therefore did a order analysis (resampled), and got the following plot: This plot shows the angular rotation plotted against amplitude.

2: Now, an FFT was done with this code:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

class FastFourierTransform:
# Amplitudes is a row vector
def __init__(self, amplitudes, t):
self.s = amplitudes
self.t = t

# Plotting in the input domain before fft
def plot_input(self):
plt.ylabel("Amplitude")
plt.plot(self.t, self.s)
plt.margins(0)
plt.show()

'''
The second half of this array of fft sequence have similar frequencies
since the frequency is the absolute value of this value.
'''
def fft_transform(self):
mean_amplitude = np.mean(self.s)
self.s = self.s - mean_amplitude # Centering around 0
fft = np.fft.fft(self.s)

# We now have the fft for every timestep in out plot.

# T is the sample frequency in the data set
T = self.t - self.t # This is true when the period between each sample in the time waveform is equal
N = self.s.size  # size of the amplitude vector
f = np.linspace(0, 1 / T, N, )  # start, stop, number of. 1 / T = frequency is the bigges freq
plt.ylabel("Amplitude")
plt.xlabel("Frequency [Hz]")
y = np.abs(fft)[:N // 2] * 1 /N

# Cutting away half of the fft frequencies.

sns.lineplot(f[:N // 2], y) # N // 2 is normalizing it
plt.margins(0)
plt.show()
time = f[:N // 2]
return fft, time

3. Result, with the normalized amplitudes plotted: Questions:

1. Does this thought process look correct?

2. Is it correct to say that the final fft plot is in the frequency domain? From this link, http://zone.ni.com/reference/en-XX/help/372416L-01/svtconcepts/svcompfftorder/, it looks like the final plot domain should be in the order domain, but I am not sure since the fft was done from the radian domain.

• Interesting. Can I please ask what your $Fs$ was?
– A_A
Jan 30 '20 at 15:43
• @A_A The resampled frequency was 1/1500. Jan 31 '20 at 8:43
• Alright, but what was the sampling frequency of the originally acquired data and what was the range of rpm?
– A_A
Jan 31 '20 at 8:46
• @A_A Oh ok, the original data was around 25kHz. The rpm range for the 10 second interval we are looking at is around 800-1200 rpm. Jan 31 '20 at 9:03

There are a few things going on here:

Does this thought process look correct?

Strictly, No.

The main reason for this is that you are not performing a resampling step that is the differentiating point between the Discrete Fourier Transform (DFT) and order analysis.

Is it correct to say that the final fft plot is in the frequency domain? From this link, ...

Yes, it is correct to say that the final FFT plot is in the frequency domain. But, the order plot would be in a quantity that has no units. It is basically a multiplication factor.

A little bit about Order Analysis and how it is different to the DFT:

OK, so you have a rotating shaft which you sample at some $$Fs$$ using an accelerometer. Let's use the figure provided, at 25kHz.

If the shaft is rotating at 1 Revolution Per Minute (RPM) then 1 revolution fits in $$60 \cdot 25000$$ samples. If you continued sampling for a full 2 minutes, you would have $$2 \cdot 60 \cdot 25000$$ samples that would also happen to describe 2 revolutions. The data for the first revolution are in indices 0-1499999, the data for the second revolution are in indices 1500000-2999999 and so on.

Now, we throttle the motor up and get the shaft at 2 RPM. Now, $$60 \cdot 25000$$ samples describe 2 revolutions! The data for the first revolution are in indices 0-749000, the data for the second revolution are in indices 750000-1499999 and so on.

More importantly, what has happened from the point of view of the data that we are interested in (the data about 1 revolution), the sampling frequency has been halved.

At 1 RPM we observed $$60 \cdot 25k$$ samples for one revolution. At 2 RPM we observed $$\frac{60 \cdot 25k}{2}$$ samples for one revolution. At 3 RPM...and so on.

At 1200 RPM, one shaft rotation takes $$\frac{1}{1200} \cdot 60 = 0.05$$ seconds and fits in 1250 samples (at 25kHz).

So what?

Imagine that there is a "bump" (for whatever reason, a bearing, uncentered load, structural resonance, whatever). At 1200 RPM, you can hear that thing banging across the floor. But at 1 RPM, the discontinuity is so "shallow" that it is spread across millions of samples, if it is even there. At low RPM, the elasticity of the materials might allow the power of this discontinuity to scale linearly. But at high RPM, the non-linearities kick in. Now, the shaft cannot bang off of its housing (it can, but by then it is too late) and you can see this because the structure will give rise to additional harmonics. So, it is not only that at high RPM things become louder, but also, new things show up, that were not observable at low RPM.

So, we have a problem here, because with increasing RPM we are trying to fit more information at increasingly less and less and less bandwidth.

Part of the answer here is to increase the sampling frequency.

But, we would still be left with no way to normalise the data for the increase in speed. Therefore, if we were to run DFT on a block of data, we would not be sure how much of a rotation or where does interesting behaviour occurs.

One approach to this would be to physically modulate the sampling frequency by the RPM.

This is what Order Analysis does, but it does it in software. What you are doing prior to DFT is resampling the signal (as a function of RPM), so that, no matter what the RPM does (going up or down), 1 shaft revolution sits in exactly the same amount of samples.

By doing this, you can observe consistent parts of the spectrum as the RPM increases (e.g. harmonics that are predictable as per the linked article) and whether interesting components start to appear at specific RPM.

A fundamental assumption of the DFT is that the signal at its input is periodic. So, although you may be passing it a voice recording of someone saying "one, two, three", what the DFT "sees" is "one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, ..." and analyses this as a whole.

This is the reason why the spectrum looks like two peaks connected with a curving ribbon on the first spectrum from the figure of the book you are posting.

What has happened in the second picture is to adjust the timing of the samples based on the RPM, so that each revolution now takes up exactly the same number of samples.

There is nothing additional hapenning to that sinusoid, except modulating its frequency in a way proportional to the RPM "signal". If you adjust for the RPM, what you see is basically just the same sinusoid.

A better way to visualise this is via the spectrogram. Here is a sweeping sinusoid: And here is its spectrogram: (Both images are from the Wikipedia article).

What the spectrogram is, is repeated applications of the FFT in short overlapping blocks throughout the signal. Therefore, its x-axis is time, its y-axis is frequency and the cell colour is proportional to the strength of that harmonic.

If you apply a spectrogram to your data, you will notice that the "run-up" (or run-down) shows up clearly as increasing (or decreasing) "lines". The reason why there might be more than one lines is because there are additional components to the rotation that also happen faster (or slower), depending on the RPM.

What we really want is to get rid of the "run-up" (or run-down) component of the DFT.

What we are interested in is the spectrum associated with 1 shaft revolution at some fixed RPM. This is equivalent to one vertical slice of the spectrogram at some time instant (which however in our case is proportional to RPM).

With order analysis, what you are doing is basically "rotating" the spectrogram, at an angle that matches this increasing line. If you did this, in this particular example, you would see a straight line. (Technically, it is not a rotation, but a "skewing", but just for illustration purposes, think of it as rotating this line so that it is parallel to the x-axis).

This is what this re-sampling (mentioned in the figure you attach) refers to.

After you do that, the x-axis of the spectrogram would be RPM, the y-axis of the spectrogram would be Order and the cell colour would be proportional to the strength of that harmonic.