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

Question

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.

Answer 1

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);
shading interp
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.

Answer 2

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

mat = loadmat('se_dsp.mat')

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

# 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[0].plot(t_sp, sp_bandpass, label='Seatpan', alpha=0.5)
axes[0].plot(t_gb, gb_bandpass, label='Gearbox', alpha=0.5)
axes[0].set_title('Acceleration near 3-4 Hz')
axes[0].set_ylabel('3-Axis Acceleration\nMagnitude (G\'s?)')
axes[0].set_ylim(-0.75, 0.75)
axes[0].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[1].plot(t_gb, gb_engine_RPM)
axes[1].set_ylim(60*(engine_Hz - 2.5*bandwidth/2), 60*(engine_Hz + 2.5*bandwidth/2))
axes[1].set_title('Gearbox Instantaneous Frequency')
axes[1].set_ylabel('Est. Engine RPM')
axes[1].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[2].plot(t_gb, gb_prop_RPM)
axes[2].set_ylim(60*(prop_Hz - 2.5*bandwidth/2), 60*(prop_Hz + 2.5*bandwidth/2))
axes[2].set_title('Gearbox Instantaneous Frequency')
axes[2].set_ylabel('Est. Prop RPM')
axes[2].set_xlabel('Time (sec)')
axes[2].set_ylim(1900, 1960)
plt.show()

Answer 3

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/