Oscillations on Hilbert envelope

Question

I am trying to extract the decay properties of an oscillating time history. The obvious method for doing this is to use a Hilbert transform. The figure shows an example of a theoretical decay and the envelope calculated using the Hilbert transform.

My problem is that the envelope has oscillations (Gibbs oscillations) due to the fact that there is a discontinuity in the sudden start of the data.

The extracted decay rate, below, which should be a constant in this theoretical case, is unacceptably corrupted by the oscillations.

I have considered filtering the envelope with a low pass filter but as the oscillations are the same order as the original frequency I will spoil the decay data. Alternatively I could create a mirror image of the data and thus remove the starting step. However, this requires finding the starting phase of the time history. Could these methods work or can you suggest another method?

Edit As suggested by Olli Niemitalo I have had another look at filtering the envelope signal. The filtered envelope is added here

The filtered envelope signal now has no oscillations but has dropped down from the peaks. This is fine for my application were I need the rate of decay. The next figure compares the decay calculated from the filtered envelope to the unfiltered envelope

This is success! Thanks to Olli.

Thanks

Answer 1

A close alternative to a Hilbert transformer is using two infinite impulse response (IIR) all-pass filters that have a phase difference of approximately 90° over the frequency range of interest. For envelope detection, the two all-pass outputs are used in place of the original signal and its Hilbert transform. Here is one such filter pair with quadrature outputs $\text{out}_0$ and $\text{out}_1$:

$$\text{mid}_0[t] = 0.4794008656\,(\text{in}[t-1] + \text{mid}_0[t-2]) - \text{in}[t-3]\\ \text{mid}_1[t] = 0.8762184935\,(\text{mid}_0[t] + \text{mid}_1[t-2]) - \text{mid}_0[t-2]\\ \text{mid}_2[t] = 0.9765975895\,(\text{mid}_1[t] + \text{mid}_2[t-2]) - \text{mid}_1[t-2]\\ \text{out}_0[t] = 0.9974992559\,(\text{mid}_2[t] + \text{out}_0[t-2]) - \text{mid}_2[t-2]$$ $$\text{mid}_3[t] = 0.1617584983\,(\text{in}[t] + \text{mid}_3[t-2]) - \text{in}[t-2]\\ \text{mid}_4[t] = 0.7330289323\,(\text{mid}_3[t] + \text{mid}_4[t-2]) - \text{mid}_3[t-2]\\ \text{mid}_5[t] = 0.9453497003\,(\text{mid}_4[t] + \text{mid}_5[t-2]) - \text{mid}_4[t-2]\\ \text{out}_1[t] = 0.9905991567\,(\text{mid}_5[t] + \text{out}_1[t-2]) - \text{mid}_5[t-2]$$

The coefficients were optimized for minimum largest error over a 20 Hz to 22.03 kHz pass band at a sampling frequency of 44.1 kHz. For a sinusoidal input the largest possible residual error in the envelope is -44 dB or ±0.6%. The following envelope is detected for a unit impulse at time 0:

The vertical axis is magnitude and the horizontal axis is time in samples. Because the all-pass filters are causal, the envelope will not be disturbed until after any transient. Reversing the exponentially decaying input signal gives quite a nice envelope all the way to the transient:

There is a bit of a time lag, which is why the input signal peaks are higher than the envelope. The time lag will be worse for a frequency close to the design band limits of the all-pass filter pair.

If your bandwidth and error requirements are different, you can optimize your own polyphase allpass filters using the HIIR library by Laurent de Soras.

About filtering the envelope: A basic property of linear time invariant (LTI) systems is that if you filter an exponential function then the filter's output is the same exponential function multiplied by a constant. That does not affect the calculated decay rate. The same applies to your truncated exponential envelope if you only use smoothing filter outputs that are not affected by the end effects, by discarding part of the output or by using an anti-causal smoothing filter. The perspective is that you would use smoothing to reduce noise in general, not to reduce end effects. But if you do use it to reduce end effects, I don't think there is an equivalence between envelope filtering and generating the envelope using an all-pass filter pair, because the first does not preserve signal power.

Answer 2

The Hilbert transform, as applied to finite (rectangularly windowed) segments of data, is only useful if your ignore startup (and/or end) transients (until the transform's impulse response drops below your desired noise floor).

You can either cut off, or smooth (or perhaps apply machine learning?) both end transients.

If your data is integer periodic in aperture width, then you can use an FFT for the transform and assume the length is infinite (no start transient).

Answer 3

The Hilbert transform is certainly the best option to extract the envelope of an oscillating signal, but as you pointed out, it may be giving oscillations when your signal deviates from the hypothesis that led to the formulation of the transform (here smoothness).

Aren't you hammering a needle with a sledgehammer?

Since you are interested in the decay rate (and that you certainly assume that it is stationary), I would recommend using a simpler method by fitting your signal with a model of it. For instance, you can write:

def model(decay, start, amp, freq, phase):
    t = np.linspace(0, N)
    signal =  (t > start) 
    signal *= np.exp(-(t-start)/decay ) 
    signal *= np.sin(2*np.pi * (freq * t -phase))
    return signal

and then use one of the zillions implementations of a least-square optimizer. In python, lmfit is my favorite. Instead of getting a result as a function of time, you'll directly get what you are looking at.

Answer 4

The solution given by @Olli Niemitalo taught me a lot (and I have not understood it completely yet). Let me try a trick without filtering (the latter may have parameters). I reduce the transient effect by mirroring the signal around its first and last "true maxima",

and computing the Hilbert envelope on the resulting extension. The same can be done around the first and last minima. Here are the obtained envelopes:

Finally, one could average the bottom and the top envelopes to further reduce oscillations (they could be post-filtered as well). This is neither very systematic, nor grounded on a solid body of work, but quite reduced in parameters. Feedback appreciated on the given Matlab code

% Design a decaying sine
nSample = 1024;
time = linspace(0,1,nSample)';
timeCorr = linspace(-1,1,2*nSample-1)';
frequency = 22;
decay = 1.34;
dataSin = sin(2*pi*frequency*time);
dataEnv = (exp(-decay*time)).^1.2;
data = dataSin.*dataEnv;

% Compute local maxima
% Uses http://www.mathworks.com/matlabcentral/fileexchange/12275-extrema-m--extrema2-m
[xmax,imax,xmin,imin] = extrema(data);
imax = sort(imax);
imin = sort(imin);

% Crop data between the first maximum and the ultimate
iMax1 = imax(1); if isequal(iMax1,1), iMax1 = imax(2);end
iMax2 = imax(end); if isequal(iMax2,nSample), iMax2 = imax(end-1);end
dataMod = data(iMax1:iMax2);
lDataMod = length(dataMod);

% Half-sample mirror symmetry at both ends
dataMod3 = [flipud(dataMod(2:end));dataMod;flipud(dataMod(1:end-1))];
timeMod3 = (iMax1-lDataMod+1:iMax2+lDataMod-1)';
figure(1);clf
hold on
plot(timeMod3,dataMod3)
plot(iMax1:iMax2,dataMod,'r')
hold off

% Hilbert enveloppe on the mirror symmetrized signal
dataHMod3 = abs(hilbert(dataMod3));
dataHTop = NaN*ones(size(data));
dataHTop(iMax1:iMax2) = dataHMod3(lDataMod:end-lDataMod+1);

%% Do the same on the lower envelope
iMin1 = imin(1); if isequal(iMin1,1), iMin1 = imin(2);end
iMin2 = imin(end); if isequal(iMin2,nSample), iMin2 = imin(end-1);end
dataMod = -data(iMin1:iMin2);
lDataMod = length(dataMod);

% Half-sample mirror symmetry at both ends
dataMod3 = [flipud(dataMod(2:end));dataMod;flipud(dataMod(1:end-1))];
timeMod3 = (iMin1-lDataMod+1:iMin2+lDataMod-1)';
figure(2);clf
hold on
plot(timeMod3,dataMod3)
plot(iMin1:iMin2,dataMod,'r')
hold off

% Hilbert enveloppe on the mirror symmetrized signal
dataHMod3 = abs(hilbert(dataMod3));
dataHBot = NaN*ones(size(data));
dataHBot(iMin1:iMin2) = dataHMod3(lDataMod:end-lDataMod+1);

figure(3);clf;hold on
plot(time,[data],'k')
plot(time,[dataEnv,-dataEnv],'kx')
plot(time,[dataHTop,-dataHBot],'r')
axis tight
legend('Data','Env. True Top','Env. True Bottom','Est. Top','Est. Bottom')
hold off