Smooth envelope of an aperiodic signal

Question

I have a signal that looks like this: Red line is the function I'd like to obtain somehow. Everything <= -60 is basically absent data (i.e. not low signal).

If I understand correctly, Hilbert transform won't work for me since the signal is aperiodic. I tried spline interpolation but it doesn't produce anything even remotely similar to what i'm looking for.

I understand this might be too vague but it's the best I could come up with :)

UPD: here's the data

Answer 1

Here is a method I came up with for a similar problem recently:

Intialize an output array of zeros with the same length as your signal.

Iterate over your list of data points. For each data point, calculate values for some smooth bell shape like shown below (orange, shown for one data point).

I used a cosine shape for this but any roughly bell-like shape could work. Specifically, I am iterating from j=-n to j=n, where n is the width of the bell, and calculating output[i+j]=input[i]*math.cos(0.5*j*(math.pi/n))

If any point on the calculated bell is higher than the current value for output at that point, overwrite your output with the higher value. i.e. output[i+j]=max(output[i+j],input[i]*math.cos(0.5*j*(math.pi/n)))

You will get something that looks like this:

You can choose to apply a linear phase lowpass filter (scipy.signal.filtfilt) to round over the sharp corners at the bottom, sacrificing some precision at the peaks.

Answer 2

UPDATED TO INCLUDE OP'S DATA AND DEMONSTRATING LOOSE AND TIGHT TRACKING

A simple way to do this is a leaky peak hold optionally followed by a low pass filter to further smooth the waveform it desired. Also see this answer which can also be applied here where an alternate peak detector was used with a Gaussian filter for a symmetrical response around the detected peaks.

The leaky peak hold is implemented by this simple conditional equation:

if $y[n]> |x[n]|$ then $y[n]=|x[n]|$, else $y[n] = \alpha y[n-1]$

The parameter $\alpha$ is the leak rate and is always <1. The closer $\alpha$ is to one, the slower the decay.

The block diagram is shown in the graphic below.

An example result is plotted below without the final low pass filter (LPF) with $\alpha =0.92$ (a lower $\alpha$ value results in faster decay and tighter tracking to short term changes):

This is subsequently smoothed by the use of a low pass filter, which can be implemented as a simple moving average or cascade of moving averages for more smoothing. The above result after subsequent filtering is shown in the plot below:

The beauty is in its simplicity and ease of adjusting the tracking as demonstrated by the python code I used for the above example:

alpha=0.92
out = [0]
for sample in wfm:
    out.append(max(np.abs(sample), alpha*out[-1]))
coeff = sig.firls(15, [0, .4, .6, 1], [1,1,0,0])
out2 = 1.06*sig.filtfilt(coeff,1,out)
plt.plot(wfm)
plt.plot(out2)

Below is the same process with "looser tracking" using $\alpha =0.995$:

Further explanation of operation is as follows:

First take the absolute value of the signal (looks like the OP has already done this).

Then test if the absolute value is greater than the output of the leaky peak hold. If so reset the output of the leaky peak hold to be this value, otherwise the output is the previous output multiplied by the leak parameter $\alpha$.

An optional low pass filter will then smooth the edges. The cutoff determines the tracking rate of the peak; the lower the cutoff the slower the allowed rate of change of the peak. A gain adjustment compensates for the averaging operation of the filter; the tighter the filter the longer the effective averaging window will be and therefore a greater gain is required. The value for the leak rate and cutoff of the filter can be determined empirically.

In my code, I used filtfilt as implementation of the low pass filter operation. This is non-causal zero-phase filter appropriate for post-processing as done here. As a zero-phase filter, the output is aligned with the input so I used it for that convenience. In actual implementation, a traditional FIR low pass filter or even a simple moving average could be used with the output inevitably delayed from the input.

Why Hilbert Won't Work

To clear up a possible misunderstanding by the OP, the reason a Hilbert approach won't work is NOT because the signal is aperiodic. It is because the Hilbert will provide an envelope of a signal at a given passband carrier frequency (with a carrier higher than the occupied signal bandwidth) which is not the case here. This is just mixing terminology in what we mean by "envelope". When we use the Hilbert to determine the envelope, we are specifically referring to the envelope as being the baseband amplitude modulation signal of a sinusoidal carrier frequency, where that frequency is higher than the highest frequency component of the baseband signal. In the OP's case, it is the baseband signal that is provided, so the result is already the "envelope" in that context.

Let me explain this by using the OP's example:

Below is the same test signal I used for the example above, with a view zoomed in on the first 2500 samples. What the OP wants is an output waveform that transitions the major peaks shown, while what is shown is the input waveform I used without any modification.

If this signal was the amplitude modulation for a higher frequency carrier, that modulated signal would appear as follows, with the result of using the absolute value of the Analytic Signal (Hilbert envelope approach) to extract the "envelope" in that context (here zoomed in on the first 600 samples of the waveform above:

Here we see the actual sinusoidal carrier with amplitude modulation (AM), with the "envelope" being the modulated amplitude waveform specifically. This is the typical use case and application for using the Hilbert Transform to extract the envelope. Other posts here may suggest that the signal needs to be narrow band and that is not the case; we only need to know the carrier frequency and that the carrier frequency is above the bandwidth of the baseband signal for this to work properly.

Bottom line, this is clearly a different envelope than what the OP is referring to, and why a leaky peak hold would be an easy solution to get what the OP really wants.

Answer 3

You may use the first step of the Empirical Mode Decomposition which interpolates the upper envelop and the lower envelop.

Basically the idea is to locate extremes and interpolate between them.
This way not all samples are effective.

In your case you're after the upper envelop.

For Python you may use the following:

• EMD Notebook - Look at the envelop.
• EMD Package - Look at the functions emd.sift.get_padded_extrema() and emd.sift.interp_envelope().

Using MATLAB, I could achieve this (On the OP's data):

You may filter it a bit go get something smoother.

This is the code for reference:

clear();
close('all');
vS = readmatrix('Data.csv');
numSamples = size(vS, 1);
minPeakProminence = 25;
minPeakDistance   = 10;
[vPeak, vPeakIdx] = findpeaks(vS, 'MinPeakProminence', minPeakProminence, 'MinPeakDistance', minPeakDistance);
[vT, vTIdx] = findpeaks(vS);
% Prevent extrapolation at the edges
if(vPeakIdx(1) ~= vTIdx(1))
    vPeakIdx = [vTIdx(1); vPeakIdx];
    vPeak    = [vT(1); vPeak];
end
if(vPeakIdx(end) ~= vTIdx(end))
    vPeakIdx = [vPeakIdx; vTIdx(end)];
    vPeak    = [vPeak; vT(end)];
end
figure();
plot(vS, 'DisplayName', 'Input Signal');
hold('all');
plot(vPeakIdx, vPeak, '*', 'DisplayName', 'Local Peak');
vC = makima(vPeakIdx, vPeak, 1:numSamples);
plot(vC, 'DisplayName', 'Interpolated Peaks');
legend();