# Oscillations on Hilbert envelope

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

• Have you tried setting the initial conditions of the Hilbert transformer, to lessen the effect of the initial discontinuity? If you're using matlab this means setting the zi input to the filter command. – Peter K. Nov 9 '15 at 12:38
• @PeterK. I am using Fourier transforms to calculate the Hilbert transform. Are you suggesting I use a filter with a 90 deg phase shift? I could see how this could work particularly if I reversed my data. However, I understand that this filter can not be completely flat. Is this correct? – Hugh Nov 9 '15 at 12:43
• Why don't you just take the local maxima of the oscillating curve to determine the rate of decay? – Matt L. Nov 9 '15 at 12:43
• @MattL. As it is sampled data there will not typically be a point at the maxima. The result is a very ragged decay plot. Also, this disregards all the points between the maxima which also contain information regarding the decay. There could be an approach involving interpolation of the data which could give better maxima but this is an alternative well away from the Hilbert approach. – Hugh Nov 9 '15 at 12:58
• @Hugh Well, the Fourier approach to calculating the Hilbert transform won't be "completely flat" either. You just don't see it because you're only looking at the sampled points. – Peter K. Nov 9 '15 at 13:17

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.

• Thanks for this. Building a causal filter and reversing the data seems to be the way to go. As my signals are narrow band there is some flexibility for my filter design. Is there any reason you used two filters ? Would one more complex filter work just as well? – Hugh Nov 10 '15 at 15:04
• Because it's a very efficient structure for the results it gives. Alternatively, with the input the first quadrature component, and a single filter as other, the filter would have to be a 90 degree phase shift, which has a non-causal anti-symmetrical impulse response. – Olli Niemitalo Nov 10 '15 at 17:28

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).

• Thank you for your comments. My data is not periodic so the last point is not a way forward. My data does taper to zero for large times so I only have to fix up one end. However, the details at the end are important so I must do this carefully. – Hugh Nov 10 '15 at 10:07

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.

• Thanks for the thought. What you are suggesting is a parametric model. Unfortunately I need a non-parametric model because my decay is non-linear. Consequently I need a good look at the envelope to detect what is happening. Even worse the start is where the non-linearity can be most significant so this is a crucial part. – Hugh Nov 14 '15 at 16:21
• I was mainly sharing a situation that recently happened to me - starting with the same analytical tools and ending up doing a simple quadratic fit with a linear model. My understanding of the Hilbert transform is that it will not be able to resolve dynamics at the finest scale. Looking forward to solutions to your problem! – meduz Nov 16 '15 at 10:15

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

• Thanks for this. I had a look at this as well as I indicated in my question. The problem is how to make a smooth mirror image. The data may not have points at a maximum and although you can find a near point this can lead to a discontinuity. I found that the best method was to interpolate the data and find a zero crossing. I then re-sampled making the zero crossing a point. Finally I did a reflection and inversion to make an extension. However, the data is not continuous in its first derivative because of the exponential decay, if considered alone, has a cusp at the join. Still this does help. – Hugh Nov 14 '15 at 22:18
• My bad, I did not see this line in your post. One advantage of mirrorring at a local extremum is that you keep the continuity in the first derivative – Laurent Duval Nov 14 '15 at 22:23
• I think the first derivative is not continuous. The exponential term will give a cusp. Your point of filtering after the join would help to smooth this out. – Hugh Nov 14 '15 at 22:28