cracking the propeller engine code

06 Mar, 2024

I took on this project to analyze the audio of an aircraft engine during takeoff (for the purpose of this post, the details of which I will keep confidential), digging into the frequencies to see what they could tell me about propeller and engine behavior. Think of it like eavesdropping on an aircraft engine, but with math and code instead of a borescope on a live engine. I’ll try to walk through everything I did, and what I found, plus some pretty cool graphs.

Loading the Audio and Breaking It Down

I started with a snippet of engine noise I wanted to dissect. I used librosa, a Python library that’s a beast for audio analysis. First step, load the audio with y, sr = librosa.load(audio_file_path, sr=None). Here, y is the audio signal (a bunch of amplitude values over time), and sr is the sample rate i.e., how many samples per second. Keeping sr=None means I let the file’s native rate (probably 44,100 Hz) take the wheel.

Next, I hit it with a Short-Time Fourier Transform (STFT) using stft = librosa.stft(y). The STFT chops the audio into short overlapping windows (default 2048 samples), applies a Fourier Transform to each, and gives you a matrix of frequency content over time. Mathematically, for a signal ( x(t) ), the STFT is:

S T F T (t, f) = \int_{- \infty}^{\infty} x (τ) w (t - τ) e^{- j 2 π f τ} d τ

where ( w(t) ) is a window function (Hann, by default), ( t ) is time, ( f ) is frequency, and ( j ) is the imaginary unit. The output stft is complex (real and imaginary parts), so I split it with stft_magnitude, stft_phase = librosa.magphase(stft). The magnitude tells me how loud each frequency is, and the phase keeps the timing info (which I’d need later).

Filtering for Propeller Frequencies

I figured the propeller’s hum would live between 50 Hz and 2000 Hz, low rumbles to higher whines. To focus there, I got the frequency bins with frequencies = librosa.fft_frequencies(sr=sr, n_fft=stft.shape[0]), which maps each row of the STFT to a frequency. The formula’s simple:

f_{k} = \frac{k \cdot s r}{n_{fft}}

where ( k ) is the bin index, ( sr ) is the sample rate, and nfft=2048n_{\text{fft}} = 2048n_{\text{fft}} = 2048 is the window size. I grabbed the bins between 50 and 2000 Hz with propeller_bins = np.where((frequencies >= 50) & (frequencies <= 2000))[0], then made a filtered STFT by zeroing out everything else: filtered_stft[propeller_bins, :] = stft[propeller_bins, :]. To turn it back into audio, I used the inverse STFT: filtered_y = librosa.istft(filtered_stft * stft_phase). This gave me a time signal with just the propeller range.

Finding Peaks and Intervals

I wanted to spot big moments in the filtered signal, so I used find_peaks from scipy.signal: peaks = find_peaks(filtered_y, height=0). This finds indices where the signal spikes above zero. Then, I calculated time intervals between peaks with intervals = np.diff(peaks) / sr, dividing by the sample rate turns sample counts into seconds. These intervals could hint at propeller rhythm, but honestly, they were a bit messy.

Spectrogram Time

Now, the fun part, visuals! I plotted a spectrogram with librosa.display.specshow(librosa.amplitude_to_db(stft_magnitude, ref=np.max), sr=sr, y_axis='log', x_axis='time'). The amplitude to dB conversion is:

dB = 20 \cdot \log_{10} (\frac{amplitude}{ref})

where ref=np.max scales it to the loudest point. This gives a heatmap of frequency vs. time, with intensity in decibels. I overlaid red lines for the 50-2000 Hz range to highlight the propeller zone.

Propeller Magnitude Peaks

I zoomed into the propeller range’s magnitude with propeller_magnitude = stft_magnitude[propeller_bins, :], then summed it across frequencies for each time step: propeller_sum_magnitude = np.sum(propeller_magnitude, axis=0). To find big shifts, I used find_peaks again, setting a threshold at mean plus two standard deviations: peaks = find_peaks(propeller_sum_magnitude, height=np.mean(propeller_sum_magnitude) + 2 * np.std(propeller_sum_magnitude)). Plotted it on the spectrogram with blue X’s for peaks.

Zeroing In on 100 Hz

Next, I picked 100 Hz, a low end frequency to analyze solo. Found its bin with specific_bin = np.argmin(np.abs(frequencies - 100)), grabbed its magnitude, and hunted peaks the same way.

Then, I graphed those peak amplitudes over time: plt.plot(time_intervals, engine_power_changes, marker='o'), assuming amplitude ties to engine power (an assumption).

Scaling Up - 50 to 500 Hz

I went bigger, analyzed 50 to 500 Hz in 50 Hz steps. For each frequency, I found the bin, extracted magnitude, spotted peaks, and plotted them. First, I tried overlaying them on a spectrogram, then graphed each frequency’s peak amplitudes over time. I did two versions: one plot with all frequencies together

and plot of each frequency:

The Math Behind the Peaks

The find_peaks function uses a simple rule: a point’s a peak if it’s higher than its neighbors and above the threshold. For a signal ( x[n] ), a peak at ( n ) satisfies:

x [n] > x [n - 1] and x [n] > x [n + 1] and x [n] > height

Here, height = mean + 2 * std is a statistical cutoff, about 95% of a normal distribution falls below it, so peaks are rare, significant blips. The amplitudes (properties['peak_heights']) are raw STFT magnitudes, which I tied to “engine power” as a rough guess.

What I Found

So, what’d I get? The spectrograms showed a busy frequency landscape, lots of action in the 50-2000 Hz range, with peaks popping up as the engine revved. The 100 Hz plot had clear spikes, suggesting moments when that frequency dominated, which maybe tied to propeller speed or throttle changes. The 50-500 Hz multi-plot was wild, peaks at different times across frequencies, hinting at unusual engine dynamics. The intervals from the filtered signal? Scattered, ranging from fractions of a second to longer gaps, but no clear pattern jumped out.

Where I’d Go Next

The peaks suggest variations in engine power for the analyzed audio range, but linking amplitude to power is a stretch without calibration data. Nevertheless, I know where the actual data comes from and the conditions under which this data has been gathered, which explains the unusual trends. The STFT worked great for spotting frequency trends, but the time resolution’s fuzzy due to window overlap. Next time, I’d maybe try wavelet transforms for better frequency-time balance. I would love to filter noise better or model propeller harmonics mathematically.