4
$\begingroup$

my signal consists of data points that represent the (vibration signal) energy measured during process cycles of a machine. Due to the hydraulic power units emptying and filling again, some periodicity is introduced into the yielding signal of energy measurements.

I am searching for a way to remove this periodicity from the data, since this effect persists over the course of all data points, however I am interested in changes in the signal that do not stem from this physical effect.

Sadly, my DSP skills are very limited, so I am hoping that someone could help me. I work with Python.

Please find attached I) an image of a data sample and II) a sample from the data.

I enter image description here

II https://pastebin.com/rgxhhmcH // Password zzxWgpFxYb

Addition: Sadly I don't have the opportunity to get the signal from the hydraulic units. The original signal was sampled with 1MHz. From these signals I calculated the energy for each process cycle. Every 1.2s one process cycle was carried out. However, there are breaks between process cycles, e. g. at ~5100 (see image).

$\endgroup$
4
  • $\begingroup$ Can you please mention the sampling frequency of the signal? (What is the time interval between two successive values in your signal?) And also, whether you have (or can get) the signal that drives the hydraulic units (?) $\endgroup$
    – A_A
    Aug 8, 2022 at 8:09
  • $\begingroup$ Hi, thank you for your response. Sadly I don't have the opportunity to get the signal from the hydraulic units. The original signal was sampled with 1MHz. From these signals I calculated the energy for each process cycle. Every 1.2s one process cycle was carried out. However, there are breaks between process cycles, e. g. at ~5100 (see image). $\endgroup$
    – derkurt
    Aug 8, 2022 at 8:31
  • $\begingroup$ Can you please confirm the following based on the Fs you provided?: Between samples 1000 and 2000 there are approximately 2 low frequency cycles. Given an Fs of 1000000 Hz, the real time interval is 0.001 seconds. If you have 2 cycles in 0.001 seconds, then your "disturbance" is approximately 2kHz. I am not sure of the application but if "hydraulics" are involved, this 2kHz sounds a bit too high. If I had to guess, I would say that your Fs is more like 10kHz. Which would make your "cycle break" at 5100 about half a second after the start. What do you think? $\endgroup$
    – A_A
    Aug 8, 2022 at 8:45
  • $\begingroup$ Sorry, I did not make the situation clear enough I think. The machine does 50 process cycles per minute. I measured a vibration signal that was sampled with 1MHz for each process cycle. Then I calculated the signal energy for each process cycle, which I have provided an image / array of. To the break is rather 1.2 s * 5100 = 6120 s after the start. $\endgroup$
    – derkurt
    Aug 8, 2022 at 9:09

1 Answer 1

6
$\begingroup$

Some time-frequency analysis with the synchrosqueezed CWT, where I guess fs=10000 but this doesn't matter except for axis labels and physical interpretation.

1. Take SSQ_CWT of original, tune wavelet for sparsity

Tune parameters until we get the most energy concentration in the fewest number of points. The signal appears to have detailed low frequency behavior, which makes CWT the suitable tool of choice as opposed to STFT.

The clear curve we see confirms this, and it's what we want to eliminate according to OP.

2. Filter the undesired

Applying SSQ inversion, which is a highly non-linear inversion that can't be achieved with FFT-based methods, we obtain

Granted, due to the simple inversion specification in this case (exclude freqs entirely), other methods may achieve similar results.

3. Take SSQ_CWT again to see what we've got, re-tune for sparsity

I also zoomed along freq. Need to change wavelet again due to the tightly-packed nature of higher frequencies in this signal (need better frequency resolution).

Zoom along time since it's hard to see:

4. Filter the undesired, round 2

The non-intense regions appear noisy - they may not be noise, but they're complicated enough that some highly nonlinear methods will be required, possibly neural nets, to make sense of them. Given that this is just machinery and not brain EEG, it's likelier noise. Zoomed along time, plus amplitude envelope:

and in SSQ:

5. Inspect amplitude via Hilbert transform

Not entirely appropriate since the single component isn't quite clean - SSQ would do better but it's a lot of work, but we should get a good approximation:

Now take its SSQ

No strong pattern, looks like noise - makes sense as the measurement is mediated by a bunch of physical processes and nothing to push the measurement intensity a specific way.

Verdict

The signal of interest is most likely an AM-FM waveform, the FM being sine-like and AM being noisy. It's considerably different from a pure sine, but depending on application, a pure sine wouldn't be a terrible approximation. A higher sampling rate would aid with further insights, as the FM bandwidth is tightly packed.

Messy code

import numpy as np
from scipy.signal import hilbert
from ssqueezepy import ssq_cwt
from ssqueezepy.visuals import plot, imshow

FS = 10000
# x = # load data

ikw = dict(abs=1, xlabel="time [samples]", ylabel="frequency [Hz]")

Tx, _, ssq_freqs, *_ = ssq_cwt(x, ('gmw', {'gamma': 1, 'beta': 10}), fs=FS)

imshow(Tx, yticks=ssq_freqs, **ikw, title="|SSQ_CWT(x)|")

x_clean = Tx[:100].sum(axis=0).real
plot(x_clean, xlabel="time [samples]",  show=1,
     title="x_clean: exclude SSQ'd freqs above %.3g Hz" % ssq_freqs[99])

#%%
Tx, _, ssq_freqs, *_ = ssq_cwt(x_clean, ('gmw', {'gamma': 3, 'beta': 30}), 
                               fs=FS, nv=256)
imshow(Tx[:500], yticks=ssq_freqs[:500], **ikw, title="|SSQ_CWT(x_clean)|",
       norm=(0, np.abs(Tx).max()/10))

imshow(Tx[:500, 512:1024], yticks=ssq_freqs[:500], **ikw, 
       xticks=np.arange(512, 1024),
       title="|SSQ_CWT(x_clean)|, zoomed", norm=(0, np.abs(Tx).max()/10))
x_cleaner = Tx[130:180].sum(axis=0).real
#%%

plot(np.arange(512, 1024), x_cleaner[512:1024], 
     title="x_cleaner: keep SSQ'd freqs between %.3g and %.3g Hz" % (
         ssq_freqs[130], ssq_freqs[180]), show=1)

#%%
Tx, _, ssq_freqs, *_ = ssq_cwt(x_cleaner, ('gmw', {'gamma': 3, 'beta': 30}), fs=FS,
                               nv=256)

imshow(Tx[:500, 512:1024], yticks=ssq_freqs[:500], **ikw, 
       xticks=np.arange(512, 1024),
       title="|SSQ_CWT(x_cleaner)|, zoomed", norm=(0, np.abs(Tx).max()/10))


xca = np.abs(hilbert(x_cleaner))
plot(xca, title="|analytic(x_cleaner)|", show=1)
#%%
Tx, _, ssq_freqs, *_ = ssq_cwt(xca, ('gmw', {'gamma': 3, 'beta': 10}), fs=FS,
                               nv=256)
imshow(Tx, **ikw, norm=(0, np.abs(Tx).max()/10), yticks=ssq_freqs,
       title="|SSQ_CWT(|analytic(x_cleaner)|)|")
$\endgroup$
5
  • 1
    $\begingroup$ Nice analysis!!! $\endgroup$
    – Peter K.
    Aug 8, 2022 at 14:37
  • 1
    $\begingroup$ Thank you very much for the insightful and detailed analysis! I have to admit that I will probably need some time to process and play a bit with your code to fully grasp what you did there. $\endgroup$
    – derkurt
    Aug 12, 2022 at 9:26
  • 1
    $\begingroup$ IMHO this answer is better written than almost all scientific papers I've seen. $\endgroup$
    – Vorac
    Jan 3, 2023 at 15:40
  • 1
    $\begingroup$ @Vorac Thanks! If you search my posts by visualization, might find more pretty stuff. $\endgroup$ Jan 6, 2023 at 13:24
  • 1
    $\begingroup$ @Vorac So I put it on steroids. $\endgroup$ Apr 13, 2023 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.