Signals and the FFT

Plot the raw output of a vibration sensor bolted to a motor and you see a fuzzy, wandering line. Hidden inside it is a clean story: the motor spins at one frequency, the gearbox adds another, and the mains supply leaks in a third. The Fast Fourier Transform (FFT) rewrites a signal as the sum of sine waves that compose it, so you can read those frequencies off directly. In Octave it is one function call, fft.

Sampling and the Nyquist limit

A microcontroller does not see a continuous voltage. It sees samples taken at a fixed rate fs, in hertz. A sine of frequency f needs more than two samples per cycle to be pinned down, so the highest frequency you can represent is the Nyquist frequency, fs / 2.

Anything faster does not vanish. It aliases: it folds back and appears at a false, lower frequency, equal to the distance to the nearest multiple of fs. A 120 Hz vibration sampled at 100 Hz shows up as a convincing 20 Hz one. This is why sensor boards put an analog low-pass filter in front of the ADC: after sampling, you cannot tell the fake from the real.

A 50 Hz plus 120 Hz test signal

We sample at 1000 Hz for exactly one second, so N = 1000 samples and the frequency spacing between FFT bins is fs / N, which is 1 Hz. Bin number k (counting from 1 in Octave) holds frequency k - 1 times that spacing.

fs = 1000;                    % sampling rate, Hz
N  = 1000;                    % number of samples (1 second)
t  = (0:N-1) / fs;            % time axis, seconds
x  = sin(2*pi*50*t) + 0.5*sin(2*pi*120*t);

amp = abs(fft(x)) / N * 2;    % single-sided amplitude
f   = (0:N-1) * fs / N;       % frequency of each bin, Hz

idx = find(amp(1:N/2) > 0.1);
printf("peak frequencies (Hz): %s\n", mat2str(f(idx)));
printf("amplitudes: %.3f %.3f\n", amp(idx));
others = amp(1:N/2);
others(idx) = 0;
printf("largest other bin: %.3f\n", max(others));

peak frequencies (Hz): [50 120]
amplitudes: 1.000 0.500
largest other bin: 0.000

Three details matter. First, fft returns complex numbers, and abs turns each into a magnitude. Second, a real signal produces a mirror image in the upper half of the spectrum, so we keep only bins up to N/2 and multiply by 2 to recover the true amplitude of each sine, which is why we get 1.000 and 0.500 as we put in. Third, the other bins are zero to three decimals because 50 and 120 hold a whole number of cycles in the record. The find call with a threshold is a crude peak picker; for noisy data the findpeaks function in Octave's signal package does it properly.

Aliasing in action

Now sample the 120 Hz sine at only 100 Hz, for 100 samples. The Nyquist limit is 50 Hz, so 120 Hz is out of range.

t2 = (0:99) / 100;
a2 = abs(fft(sin(2*pi*120*t2))) / 100 * 2;
[m, k] = max(a2(1:50));
printf("apparent frequency: %d Hz, amplitude %.3f\n", k-1, m);

apparent frequency: 20 Hz, amplitude 1.000

The spectrum shows a perfect 20 Hz tone with the full amplitude. Nothing in the data warns you.

Windowing and leakage

The FFT silently assumes your record repeats forever. If it does not contain a whole number of cycles, the end does not meet the start, the repeat has a jump, and energy leaks into every bin. Real vibration is almost never an exact bin frequency, so try 50.3 Hz, with and without a Hann window that tapers the record to zero at both ends. Dividing by 0.5, the window's average, restores the amplitude scale. Continue in the same session:

x = sin(2*pi*50.3*t);                    % between two bins
w = 0.5 - 0.5*cos(2*pi*(0:N-1)/N);       % Hann window
names = {"rectangular", "hann"};
sigs  = {x, x .* w};
gains = [1 0.5];
for n = 1:2
  a = abs(fft(sigs{n})) / N * 2 / gains(n);
  [m, k] = max(a(1:N/2));
  printf("%s: peak %d Hz, amplitude %.3f, far leakage %.4f\n", ...
         names{n}, k-1, m, max(a(201:N/2)));
end

rectangular: peak 50 Hz, amplitude 0.858, far leakage 0.0021
hann: peak 50 Hz, amplitude 0.943, far leakage 0.0000

The unwindowed peak reads 0.858 instead of 1.0, and bins 150 Hz or more away still carry 0.0021 of leakage. With the Hann window the peak is 0.943 and the far leakage is under 0.00005. The price is a slightly wider peak. This trade is worth it when a small fault signature sits next to a large motor tone: leakage from the big peak would otherwise bury it.

The same maths in NumPy

numpy_version.py

ResetRun

import numpy as np fs = 1000 # sampling rate, Hz N = 1000 # samples (1 second) t = np.arange(N) / fs x = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t) amp = np.abs(np.fft.fft(x)) / N * 2 # single-sided amplitude f = np.arange(N) * fs / N # frequency of each bin half = amp[:N // 2] idx = np.where(half > 0.1)[0] print("peak frequencies (Hz): [%s]" % " ".join("%g" % v for v in f[idx])) print("amplitudes: %.3f %.3f" % tuple(half[idx])) others = half.copy() others[idx] = 0 print("largest other bin: %.3f" % others.max()) # Aliasing: a 120 Hz sine sampled at only 100 Hz t2 = np.arange(100) / 100 a2 = np.abs(np.fft.fft(np.sin(2 * np.pi * 120 * t2))) / 100 * 2 k = int(np.argmax(a2[:50])) print("apparent frequency: %d Hz, amplitude %.3f" % (k, a2[k]))

Notice that NumPy bins start at index 0, so bin 50 is exactly 50 Hz and no minus one is needed. The same listing with a window:

numpy_windowing.py

ResetRun

import numpy as np fs, N = 1000, 1000 t = np.arange(N) / fs x = np.sin(2 * np.pi * 50.3 * t) # between two bins w = 0.5 - 0.5 * np.cos(2 * np.pi * np.arange(N) / N) # Hann window for name, y, gain in [("rectangular", x, 1.0), ("hann", x * w, 0.5)]: a = np.abs(np.fft.fft(y)) / N * 2 / gain k = int(np.argmax(a[:N // 2])) far = a[200:N // 2].max() print("%s: peak %d Hz, amplitude %.3f, far leakage %.4f" % (name, k, a[k], far))

Change 50.3 to 50.0 and both rows read 1.000: leakage only appears when the frequency does not sit on a bin. In the aliasing lines of the first runner, change 120 to 130 and the apparent frequency moves to 30 Hz, the distance to the nearest multiple of 100.