# What algorithm and approach could isolate low-frequency vibration in an engine-driven, double gear-reduced propeller system?

I have a drivetrain setup where I have a ~3Hz vibration which I can feel quite strongly in the seat of my pants, but can't seem to see in accelerometer data.

Below is a diagram of the system. There is a 4-cycle 4-cylinder Rotax 914 engine which drives a 2.45:1 gearbox which in turn drives a driveshaft which in turn drives another 1.1:1 belt drive which finally drives the three-bladed propeller: The engine spins around 5600rpm at max speed, and at this point the prop is around 2100RPM (~2.6:1 overall reduction). At a range of engine speeds from 4200-5400rpm I feel a sharp 3-4Hz impulse through the seat of my pants. Via a series of tests we have conclusively demonstrated that the impulse is tied to engine RPM, and nothing else. However, there is absolutely nothing which spins at 3-4Hz, or anywhere close to that. And the vibration completely smooths out below 4200rpm and above 5400 rpm I can no longer feel it (although there's a lot of broadband vibration at that rpm).

This impulse's band-limited nature points toward some kind of poorly damped constructive resonance between the various spinning elements. But before concluding that, I'd like to make sure that the data supports the theory. I have hooked up an accelerometer near the seatpan, another accelerometer on the gearbox, and a tachometer to measure engine position (which can be used to get RPM). I'd like to see if the signal is tied to some order of rpm.

I have done some basic spectrometer work (with Matlab), and I cannot see any discernible low-speed signal. In the accelerometer data I can make out the propeller RPM, the engine RPM, and (I think) even the driveshaft RPM. And, yet, this signal I can so easily feel in my body I cannot see in my data. Investigating the lower frequencies shows just noise.

What is a recommended approach to extracting this low speed signal in a noisy environment?

P.S. I am using two GY-521 accelerometers, which are the knock-off versions of the MPU-6050. I measure them at 1000Hz, using interrupts so I am sure that the data is read in a timely manner.

# Update

Thanks to the people who answered, I learned that it's a beat frequency I'm feeling. The equation for a beat frequency is simply $$f_{beat} = |f_2 - f_1|$$. So since the driveshaft is 10% faster than the propeller, if the propeller spins at 2000rpm, then |(2000/60*1.1 - 2000/60)| = 3.3Hz. That's precisely what I feel.

Interestingly, that means there are two other beat frequencies as well, between the prop/engine and the driveshaft/engine. These frequencies are seen in the analysis presented in https://dsp.stackexchange.com/a/88983. My guess is that the 7.5Hz and 11Hz peak have to do with beat frequencies between cylinder-level vibration events, which occur at half engine speed, and the propeller and driveshaft, respectively.

• Are you in the seat when performing the test? Resonance frequency is inversely related to mass and could be localized to a response of the chair.
– Ash
Aug 1 at 1:03
• @Ash that's an excellent point. I am in the seat when I feel it. Someone sitting in the other seat can also feel it at the same time, so it's not localized to just the one seat pan. It's also worth mentioning that this is a very reclined seat, and not set up like a chair: th.bing.com/th/id/…. The bottom of the seat pan is an integral part of the airframe, so it is not shock-absorbed or otherwise isolated. Aug 1 at 14:38
• Should the takeaway be that I would stand a better chance of detecting the low-frequency vibration if I were to attach a second accelerometer to the seat pan? Aug 1 at 14:39
• I'm just verifying that when you take measurements with the accelerometer, you are recreating the conditions in which you experience the vibration (e.g., by sitting in the seat you may apply a specific loading to the airframe). If you feel the vibration, your sensor almost certainly does as well. Low frequencies are somewhat agnostic to mounting material (wax vs. epoxy). It could be beneficial to post the 3 channel (one for each accelerometer axis) to let us else play around with it. Have you verified that the sensor works as expected outside of your test setup?
– Ash
Aug 1 at 15:18
• @Ash I have some data, but based on what you're saying there's a good chance that my measurement points where wrong. I will go for another flight with one sensor in my seat pan and the other one on the gearbox. Aug 1 at 18:15

I've done a FFT analysis of the seatpan vertical acceleration data: After zooming in on the lower frequencies it seems that higher amplitudes occur closer to 0 Hz but not at 4 Hz: However the two increased frequencies at ~34 Hz and ~39 Hz together with the frequency at ~46 Hz should lead to some beating phenomena at low frequencies due to modulation which might be perceived as low frequency vibration.

Calculating a modulation spectrum vs time shows this: Averaging of the modulation spectra gives: Further reading about the metrics used: Modulation Analysis Vibration Analysis of Rotating Machinery (Matlab)

Edit

As I have not done the analysis with matlab I can't simply post the code. But I can briefly say what would be needed. For the FFT colorplot, given y is the data and t the time vector:

FFTlen=1; % FFT length in seconds
[~,b]=min(abs(t-FFTlen));
if rem(b, 2) ~= 0
b=b-1;
end
FrameSize = b;
Fs = 1/(t(2)-t(1));

hannEstimator = dsp.SpectrumEstimator('PowerUnits',"dBm",...
'Window',"Hann",'FrequencyRange',"onesided",...
'SpectralAverages',numAvgs,'SampleRate',Fs);

step=round(0.1*FrameSize);
FFT=[];
i=0;

while (i*step+FrameSize)<=length(y)
x=y((i*step+1):(i*step+FrameSize));
Pse_hann = hannEstimator(x);
tf(i+1)=t(i*step+1);
FFT(:,i+1)=Pse_hann;
i=i+1
end

pFFT=pcolor(tf,[0:((Fs/2)/(size(FFT,1)-1)):(Fs/2)],FFT);
xlabel('Seconds')
ylabel('Frequency (Hz)')
title('FFT')
colormap(jet(256));
colorbar
zMin=min(min(FFT));
zMax=max(max(FFT));
caxis([zMax-40, zMax]);


For the modulation analysis you need to similiarly calculate a FFT at fixed timestaps, but then from the envelope of the time signal. For calculating the envelope first you probably can use the Matlab envelope() function. Matlab envelope function

Lastly for the analysis of your beating frequencies I made this table: It's very interesting that the 1/2 engine order seems to be the cause for the 2 higher beating frequencies. I think you are right that it has to do with between cylinder variations. I once had an issue with a high half order at a vehicle engine and changing the injectors to ones with lower tolerances helped significantly. It's also interesting that you don't feel the ~8 and ~11Hz but the ~3.5Hz. So maybe it would be possible to improve that ~3.5Hz beating issue by changing the distribution of the transmission ratio from 2.45&1.1 to something like 2.25&1.2. That way the overall ratio of 2.45*1.1=~2.7 would be kept (1.2*2.25=~2.7), but the resulting beating frequency would increase to |(2000/60*1.2 - 2000/60)| = 6.7Hz which might be enough to not be felt that much anymore.

• Fascinating! The 3.3Hz signal has got to be what I was feeling. What's so interesting here is that there're stronger peaks at 7.5 and 11Hz. I don't feel those, but that might be because they're too fast to individually sense. Aug 8 at 13:49
• I'll read those links, but for completeness would you mind posting the matlab code you used? Aug 8 at 13:49
• With some further thought, and looking at my graph from engineering.stackexchange.com/questions/55917/…, I think the 7.5Hz and 11Hz signals are coming from the half-speed of the engine, which is the individual cylinder firing speed. (There's possibly a cylinder which is not firing as well as it should.) When this beat frequency is calculated for the propeller and driveshaft, it winds up in the ~8 and ~11Hz, respectively. Aug 8 at 17:34

I performed some analysis on the data you provided. The first plot is the 3-axis acceleration magnitude for both the seatpan and gearbox, sampled at 1 kHz. Both accelerations follow each other fairly consistently and, assuming the units are accurate, there are some substantial oscillations at low frequencies - so you likely aren't crazy! Note that the large oscillations at the beginning are a transient artifact of the IIR butterworth filter and should just be ignored. The stable region after about 12 seconds looks more trustworthy, so focus on that for the three plots.

The second and third plots are the instantaneous frequency estimates using the Hilbert transform about the bandpassed engine and propeller rotational frequencies, respectively. It looks like there might be some resonances closer to 5175 RPM that are leading to your seatpan vibration. It's difficult to make many conclusions from this data about where the vibration is propagating from, but hopefully this helps you track it down.

I've included the code down below if you want to play around with it. Though it's in Python, the MATLAB code would look very similar. from scipy.io import loadmat
from matplotlib import pyplot as plt
from scipy.signal import butter, sosfilt, hilbert
import numpy as np

# Load the gearbox and seatpan acceleration data
gb = mat['accel_gearbox']
sp = mat['accel_seatpan']

# Sample rate (Hz)
fs = 1000

# Calculate magnitude of 3 axes of acceleration
gb_mag = np.sqrt(np.sum(gb[:, 1:]**2, axis=1))
sp_mag = np.sqrt(np.sum(sp[:, 1:]**2, axis=1))

# Create two plots that have a linked x-axis
fig, axes = plt.subplots(3, 1, sharex='all')

# Define a nominal engine frequency to bandpass filter
engine_RPM = 5100

# Convert RPM to Hz for engine and prop
engine_Hz = engine_RPM / 60
prop_Hz = engine_Hz / 2.6

#### Plot magnitude of 3-4 Hz Acceleration ####
# Perform a 10th order butterworth bandpass about 2-5 Hz
sos = butter(10, (2, 5), 'bandpass', fs=fs, output='sos')
gb_bandpass = sosfilt(sos, gb_mag)
sp_bandpass = sosfilt(sos, sp_mag)

# Get time vectors for each since they are of different lengths
t_gb = np.arange(0, len(gb_bandpass)) / fs
t_sp = np.arange(0, len(sp_bandpass)) / fs

# Plot it
axes.plot(t_sp, sp_bandpass, label='Seatpan', alpha=0.5)
axes.plot(t_gb, gb_bandpass, label='Gearbox', alpha=0.5)
axes.set_title('Acceleration near 3-4 Hz')
axes.set_ylabel('3-Axis Acceleration\nMagnitude (G\'s?)')
axes.set_ylim(-0.75, 0.75)
axes.legend()

#### Plot instantaneous frequency of engine ####
# Filter acceleration +/- 2.5 Hz about 5100 RPM (85 Hz) using 10th order butterworth bandpass
bandwidth = 5
sos = butter(10, (engine_Hz - bandwidth/2, engine_Hz + bandwidth/2), 'bandpass', fs=fs, output='sos')
gb_engine_bandpass = sosfilt(sos, gb_mag)

# Get analytical signal of filtered data using hilbert transform, perform window to attempt to manage edge effects
gb_engine_analytical = hilbert(gb_engine_bandpass*np.hanning(len(gb_engine_bandpass)))

# Instantaneous frequency is the derivative of the unwrapped phase of the analytical signal
gb_engine_RPM = np.diff(np.unwrap(np.angle(gb_engine_analytical))/(2*np.pi)*fs*60, prepend=0)

# Plot it
axes.plot(t_gb, gb_engine_RPM)
axes.set_ylim(60*(engine_Hz - 2.5*bandwidth/2), 60*(engine_Hz + 2.5*bandwidth/2))
axes.set_title('Gearbox Instantaneous Frequency')
axes.set_ylabel('Est. Engine RPM')
axes.set_ylim(5050, 5275)

#### Plot instantaneous frequency of prop, should follow engine frequency relatively closely ####
# Filter acceleration +/- 1 Hz about geared-down engine frequency using 10th order butterworth bandpass
bandwidth = 2
sos = butter(10, (prop_Hz - bandwidth/2, prop_Hz + bandwidth/2), 'bandpass', fs=fs, output='sos')
gb_prop_bandpass = sosfilt(sos, gb_mag)

# Get analytical signal of filtered data using hilbert transform, perform window to attempt to manage edge effects
gb_prop_analytical = hilbert(gb_prop_bandpass*np.hanning(len(gb_prop_bandpass)))

# Instantaneous frequency is the derivative of the unwrapped phase of the analytical signal
gb_prop_RPM = np.diff(np.unwrap(np.angle(gb_prop_analytical))/(2*np.pi)*fs*60, prepend=0)

# Plot it
axes.plot(t_gb, gb_prop_RPM)
axes.set_ylim(60*(prop_Hz - 2.5*bandwidth/2), 60*(prop_Hz + 2.5*bandwidth/2))
axes.set_title('Gearbox Instantaneous Frequency')
axes.set_ylabel('Est. Prop RPM')
axes.set_xlabel('Time (sec)')
axes.set_ylim(1900, 1960)
plt.show()

• Fantastic! I'll give it a good look. I also noticed that I felt the pulse the most strongly in my seatback, so I'll maybe move the sensor there and rerun your code. Aug 8 at 3:24

For your case of low frequency exploration, I would suggest to use autocorrelation (better than FFT for low-frequencies).

Some python examples can be found here:

https://www.geeksforgeeks.org/how-to-calculate-autocorrelation-in-python/

• That's a neat idea. As a heavy Matlab user for control systems, I have a certain familiar with cross-correlation methods, although I must admit I don't see how they would apply here. Could you expand the answer with some strategies to apply to the real-world noisy data? Aug 1 at 14:42
• I'm also interested in the comment about "low frequencies". What defines "low" in this context? Is it something like a rule of thumb that if you have a strong signal at X frequency, then, e.g., two orders of magnitude slower would be hard to see with an FFT? Aug 1 at 14:44
• When you are working with low frequencies, i.e. working near the origin in the FFT frequency domain, it can be difficult to exactly identify peaks due to a lack of resolution. For these cases, I usually rely on autocorrelation. In Matlab you can use the xcorr function. Aug 1 at 15:07
• I use xcorr for determining phase delays, but haven't used it to identify peaks in noisy data. Would you be willing to flesh out the answer to demonstrate this technique? I can play around with xcorr, but I don't have a sense as to TLAR and TLAW for real-world vibration data. Aug 1 at 15:12
• I could add some code but it would be exactly as the examples since I don't know your signals. XCORR multiplies the signal by itslef with a given lag corresponding to Ts=1/Fs. When you calculate autorcorrelation you will usually find a very high peak at 0 (lag). The next peak will be the first "hidden" periodocity. To find the peaks you can use FINDPEAKS. Aug 1 at 22:46