Time synchronous averaging - MATLAB

Question

I have a time domain signal that corresponds to a vibration signal of a machine. I have a second signal that corresponds to a tachometer signal (There is a pulse every one revolution of the shaft). I'd like to use time synchronous averaging under MATLAB. I'd like to take all blocks of time from the time domain signal whose durations is the duration of one revolution of the signal and average them out. The problem I have is that the speed of the shaft is slightly fluctuating. I've set it to 900RPM but it goes from 895RPM to 904RPM during the acquisition and so the duration of one revolution fluctuates as well (some miliseconds). My question is would averaging blocks of signals of different durations be a problem? And if so, how should I deal with this?

Answer 1

The answer to your question "is it a problem" is really up to you to determine. Sometimes the simple things work well enough.

There is a technique called Dynamic Time Warping:

https://en.wikipedia.org/wiki/Dynamic_time_warping

and Matlab has a function dtw

https://www.mathworks.com/help/signal/ref/dtw.html

which can use to figure out if it is a problem.

I suggest you pick a reference interval, that would be the closest to the average shaft interval and use that to warp the rest against.

Answer 2

In this type of problem the common approach is to resample the data from the time domain to the angle domain. The steps involved are:

1. Use the tachometer signal to obtain a speed signal. This is the inverse of the difference in time for each trigger point.
2. Depending upon the variability of the tachometer signal least squares cubic spline fitting is often used. Matlab spline library will work as will matlab central libraries. I use the fastBspline library.
3. Integrate the speed signal to obtain shaft angle position.
4. Determine the equal angle increments that you want to use. Keep the Nyquist theorum in mind, in other words there should be roughly the same number of points in the angle domain as in the time domain.
5. Switch the (time, angle) data to (angle, time).
6. Resample the (angle, time) data to obtain uniform angle increments.

Once this is done, do all the same processing you would have done, but now you are in the angle order domain.

Answer 3

Why not using the tsa in Matlab ?

Answer 4

You sample your data at fixed sampling rate in time?

However, the time delay between 2 revolutions is not constant as you mentionned your RPM is not constant.

If you split the data in blocks that correspond to one revolution, your blocks won't have the same length in time.

Answer 5

Hopefully you've found a solution or implemented your own -- if so, please share it! I'd love to see a python implementation of TSA!

Approach

So, I think the best overall approach here is to first Order Track your signal (since you have the encoder signal as well) which will leave you with a resampled signal in the angular/rotation domain. The resampling step involves a bit of interpolation as you need samples equidistant in the angular domain and not in the time domain. Because of that however, you can then split your signal in any integer number of segments and average them together, getting you the results of TSA.

Incidentally, taking the FFT of that gets you to the Order domain. Note that the frequencies of interest will now be orders, i.e. there is no need to multiply by the speed/frequency of the shaft (if you are interested in things like gear meshing frequency for example).

Note that there is no need to convert it to a rpm signal, as you would need to convert back to revolutions for the TSA anyway.

A few thoughts for implementation

Pulse detection and getting the rotations vector/array.

Now, given that your encoder pulses behave nicely from the graph you provided, it's should be easy to detect the pulses, perhaps with something like this:

import numpy as np

pulse_detection_high = 0.99 # close to 1 from looking at your graph
pulse_detection_low = 0.8 # tune this and above values to avoid false positives


# create mask of pulse locations
pulse_mask = [True if (
    value > pulse_detect_limit_high and
    encoder[index - 1] < pulse_detect_limit_low
    )
    else False
    for index, value in enumerate(encoder)
]

# convert the above to proper numpy array
pulse_indices_array = np.array(pulse_indices)

# take advantage of of True=1 and False=0
# to get a running count of pulses
pulse_cum_sum = np.cumsum(pulse_indices_array[pulse_indices])

# Convert pulses to rotations
rotations = pulse_cum_sum / cycles_per_revolution

# Corresponding pulse arrival times for rotations
rotations_time = time[pulse_indices]
# assume time above is the array accompanying your signal

# Equidistant values in rotation domain
rotations_for_interpolation = np.linspace(
    start=rotations[0],
    stop=rotations[-1],
    num=rotations.size
)

Next steps

Having calculated rotations and rotations_time above, you can link time domain to rotation domain and make any interpolations as necessary. You could create an interpolator from scipy and then feed it the rotations_for_interpolation to get the times out for when equidistant rotations have taken place and then feed that to an interpolator between signal and time to get signal values for the times you determined previously. You'll have to decide what to do with the data points before the first pulse and after the last one.

Congratulations, you have essentially order tracked your signal and now you should have a pair of numpy arrays of signal_order_tracked and rotations_for interpolation. You can now split this up into segments of the same length (that corresponds to the same angles now) and simply average over them to get TSA results.

Taking the DFT/FFT of that should give you nice indications of any deterministic signals (e.g. gears), just remember that the length of the segment will determine your order resolution (just like for the frequency domain).

Hope this helps a bit!