# Envelope detection in MATLAB

I have a sine wave of peak value 1000, frequency 1kHz sampled at 16kHz. I need to find the envelope of this signal using hilbert transform in MATLAB.I have used the inbuilt function abs(Hilbert(input signal)).This is taken as ideal case. Now I tried to find out the envelope by taking a Hilbert filter of order 63.The real part of the analytical signal is obtained by applying a delay of 31 samples and the imaginary part is obtained by passing through a Hilbert filter of order 63.Then the envelope is found out by taking the absolute value of the analytical signal. The result so obtained is compared with the ideal case output obtained earlier. There I am getting an oscillating waveform. Is it because I am using a Hilbert filter of order 63.

Is it because I am using a Hilbert filter of order 63 ?

Yes. The Hilbert Transform is difficult to implement in practice. Its impulse response is infinitely long and infinitely non-causal. That means you have to truncate, and that truncation generates the type of artifacts you are observing.

There are many other ways to determine the envelope of a signal. The best one depends on the type of signal and the requirements of you application, but it many cases there are better options than the Hilbert Transform.

The Hilbert is not at all difficult to implement in practice for narrow band signals (we just need to introduce a perfect quadrature phase shift and maintain the amplitude balance-and with a single tone as the OP has done, this can be accomplished perfectly down to our numerical precision). But for wideband signals the Hilbert is indeed challenging.

A dominant reason for the ripple in this case is NOT because of the length of the filter but because the delay is not properly matched (in particular there is an additional half sample delay introduced in the assumed delay of the Hilbert filter). The truncation will cause ripple in frequency, however not in time- it is quadrature and amplitude mismatch that will cause ripples in time, and the delay leads to quadrature or phase error. For this simple case of a narrow band tone in particular, amplitude and phase matching can be accomplished with even shorter length filters.

The Hilbert as computed is half the length of the filter, and since an odd length filter is used, there will be a half sample residual delay. For cases when we can sample at twice the rate, the proper half sample delay correction can be introduced by delaying by an integer number of samples at the higher rate and then after decimating by two will result will be a perfect half sample delay at all frequencies.

To understand how the half sample delay results in the ripple (and similarly if the amplitude was not balanced), note that the half sample delay will be a phase shift at any given frequency, and from the equations and graphics below we can see how ripple in the time domain would be introduced:

First consider the analytic signal (which is the sum of the signal plus the Hilbert transform):

$$x_a(t) = x(t) + j\hat x(t)$$

For a sinusoid this is given as:

$$x_a(t) = \cos(\omega t) + j\sin (\omega t)$$

And from this case with Euler's equation and the graphic below, we see how this reduces to $$x_a(t) = e^{j\omega t}$$; a single rotating phasor on the complex plane whose magnitude (envelope) is always 1.

Now if we were to delay $$x(t)$$ relative to $$\hat x(t)$$, this would be a fixed phase at one particular frequency, and the result would be two phasors in counter rotation, and for small angles it would appear as in the graphic below with a larger phasor added to a smaller phasor in counter rotation (consistent with a half sample offset of a 1 KHz sinusoid at a sampling rate of 16 KHz, consider the normalized radian frequency is $$\omega = 2\pi(1e3)/16e3 = 0.392$$ radians per sample, so a half sample delay in the OP's case is 0.196 radians).

The formula in the graphic is derived as follows using $$y(t)$$ as our faulty analytic signal with phase introduced:

$$y(t) = \cos(\omega t + \phi) + j \sin (\omega t + \phi)$$

Expanding into exponentials using Euler's formula:

$$= \frac{e^{j(\omega t + \phi)}+e^{-j(\omega t + \phi)}}{2} + j\frac{e^{j\omega t}-e^{-j\omega t}}{2j}$$

$$= \frac{e^{j\phi/2}e^{j(\omega t + \phi/2)}+e^{-j\phi/2}e^{-j(\omega t + \phi/2)}}{2} + \frac{e^{-j\phi/2}e^{j(\omega t + \phi/2)}-e^{j\phi/2}e^{-j(\omega t+\phi/2)}}{2}$$

$$= \frac{e^{j\phi/2}+e^{-j\phi/2}}{2} e^{j(\omega t + \phi/2)} + \frac{e^{-j\phi/2}-e^{j\phi/2}}{2} e^{-j(\omega t + \phi/2)}$$

$$= \cos(\phi/2)e^{j(\omega t + \phi/2)}-j\sin(\phi/2)e^{-j(\omega t + \phi/2)}$$

If we were to observe the smaller phasor relative to the larger phasor we immediately see how the addition of the two phasors results in envelope ripple in time! This can be compared to the actual results to see if delay is the only factor (or if there was an amplitude error as well).

Note: An FFT on the data would reveal that this is the case: the analytic signal has positive frequency components only (meaning all bins above N/2 should be zero). Quadrature imbalance in phase and amplitude causes "side-band leakage" which means the negative frequencies (those bins above N/2) will be non-zero. This is indeed how the Matlab hilbert function works (which returns the analytic signal, not actually the Hilbert transform alone but the imaginary component will be the Hilbert): the hilbert function under the hood takes an FFT of the data and zeros the negative frequencies and returns the inverse FFT of that.