First we import some packages.

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
sns.set_style("white")
sns.set_style("ticks")

Next, we define functions to create envelope and carrier signals.

def create_envelope(num_samples,fL_hz=100,N=1000,beta=3):
    # we will create the envelope by low pass filtering white noise
    mean = 0
    std = 1
    taps = signal.firwin(N, fL_hz/nyq_rate, window=('kaiser', beta))
    x = np.random.normal(mean, std, size=num_samples)
    x = signal.filtfilt(taps, 1,x)
    x = (x-np.min(x))
    x = x/np.max(x)
    return x
def create_carrier_harmonic(fs=8e3,fc=200,num_samples=1000,ncomps=1):
    # this will create a harmonic carrier with ncomps harmonics
    x = []
    for i in range(ncomps):
        x.append(np.sin(2*np.pi*(i+1)*fc*np.arange(0,num_samples,1)/fs))
    x = sum(x)
    x = x/np.max(np.abs(x))
    return x
def create_carrier_noise(mean=0,std=1,num_samples=1000):
    # this will create a white noise carrier
    x = np.random.normal(mean, std, size=num_samples)
    x = x/np.abs(x)
    return x

Amplitude modulated tone

Let's start with an amplitude modulated tone and estimate the envelope.

# init filter and signal params
fs = 8e3
dur = 5 # signal duration
num_samples = int(dur*fs)
t = np.arange(0,num_samples,1)/fs
nyq_rate = fs / 2.0
width = 5.0/nyq_rate # 5 Hz filter transition width.
ripple_db = 60.0 # stop band attenuation
N, beta = signal.kaiserord(ripple_db, width)
# create envelope
fL_hz = 50
envelope = create_envelope(num_samples=num_samples,fL_hz=fL_hz,N=N,beta=beta) 
# create carrier
carrier = create_carrier_harmonic(fs=fs,fc=100,num_samples=num_samples,ncomps=1)
# create am-fm signal
x = np.multiply(envelope,carrier)
# estimate analytic signal
ax = signal.hilbert(x)
envelope_hat = np.abs(ax)
# plt.plot(x)
fig = plt.subplots(figsize=(8,4))
ax = plt.subplot(1,1,1)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope,color='blue',label='ENVP. TRUE')
ax.plot(t,envelope_hat,color='red',label='ENVP. EST.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.xlim(0,0.5)
plt.ylim(-1,1.75)
plt.show()

Amplitude modulated harmonics

Let's make the carrier a sum of 2 harmonics. The carrier is not a tone but a sum of two tones.

# init filter and signal params
fs = 8e3
dur = 5 # signal duration
num_samples = int(dur*fs)
t = np.arange(0,num_samples,1)/fs
nyq_rate = fs / 2.0
width = 5.0/nyq_rate # 5 Hz filter transition width.
ripple_db = 60.0 # stop band attenuation
N, beta = signal.kaiserord(ripple_db, width)
# create envelope
fL_hz = 50
envelope = create_envelope(num_samples=num_samples,fL_hz=fL_hz,N=N,beta=beta) 
# create carrier
carrier = create_carrier_harmonic(fs=fs,fc=100,num_samples=num_samples,ncomps=2)
# create am-fm signal
x = np.multiply(envelope,carrier)
# estimate analytic signal
ax = signal.hilbert(x)
envelope_hat = np.abs(ax)
# plt.plot(x)
fig = plt.subplots(figsize=(8,4))
ax = plt.subplot(1,1,1)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope,color='blue',label='ENVP. TRUE')
ax.plot(t,envelope_hat,color='red',label='ENVP. EST.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.xlim(0,0.5)
plt.ylim(-.5,1)
plt.show()

The envelope estimated is a poor fit to the true envelope. Can we improve this? One option is to apply a lowpass filter to the estimated envelope. Lets see how does it do.

# init filter and signal params
fs = 8e3
dur = 5 # signal duration
num_samples = int(dur*fs)
t = np.arange(0,num_samples,1)/fs
nyq_rate = fs / 2.0
width = 5.0/nyq_rate # 5 Hz filter transition width.
ripple_db = 60.0 # stop band attenuation
N, beta = signal.kaiserord(ripple_db, width)
# create envelope
fL_hz = 50
envelope = create_envelope(num_samples=num_samples,fL_hz=fL_hz,N=N,beta=beta) 
# create carrier
carrier = create_carrier_harmonic(fs=fs,fc=100,num_samples=num_samples,ncomps=2)
# create am-fm signal
x = np.multiply(envelope,carrier)
# estimate analytic signal
ax = signal.hilbert(x)
envelope_hat = np.abs(ax)
taps = signal.firwin(N, fL_hz/nyq_rate, window=('kaiser', beta))
envelope_hat_filt = signal.filtfilt(taps, 1,envelope_hat)

fig = plt.subplots(figsize=(16,4))
ax = plt.subplot(1,2,1)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope,color='blue',label='ENVP. TRUE')
ax.plot(t,envelope_hat,color='red',label='ENVP. EST.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.xlim(0,0.5)
plt.ylim(-.5,1)


ax = plt.subplot(1,2,2)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope,color='blue',label='ENVP. TRUE')
ax.plot(t,envelope_hat_filt,color='red',label='ENVP. EST. FILT.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.xlim(0,0.5)
plt.ylim(-.5,1)


plt.show()

Amplitude modulated white noise

We do better with filter but still not perfect. Now, lets make the carrier a broadband white noise instead of tones.

# init filter and signal params
fs = 8e3
dur = 5 # signal duration
num_samples = int(dur*fs)
t = np.arange(0,num_samples,1)/fs
nyq_rate = fs / 2.0
width = 5.0/nyq_rate # 5 Hz filter transition width.
ripple_db = 60.0 # stop band attenuation
N, beta = signal.kaiserord(ripple_db, width)
# create envelope
fL_hz = 50
envelope = create_envelope(num_samples=num_samples,fL_hz=fL_hz,N=N,beta=beta) 
# create carrier
carrier = create_carrier_noise(mean=0,std=1,num_samples=num_samples)
# create am-fm signal
x = np.multiply(envelope,carrier)
# estimate analytic signal
ax = signal.hilbert(x)
envelope_hat = np.abs(ax)
taps = signal.firwin(N, fL_hz/nyq_rate, window=('kaiser', beta))
envelope_hat_filt = signal.filtfilt(taps, 1,envelope_hat)

fig = plt.subplots(figsize=(16,4))
ax = plt.subplot(1,2,1)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope,color='blue',label='ENVP. TRUE')
ax.plot(t,envelope_hat,color='red',label='ENVP. EST.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.xlim(0,0.5)
plt.ylim(-2,2)


ax = plt.subplot(1,2,2)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope,color='blue',label='ENVP. TRUE')
ax.plot(t,envelope_hat_filt,color='red',label='ENVP. EST. FILT.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
plt.xlim(0,0.5)
plt.ylim(-2,2)


plt.show()

Again, we see that filtered envelope is a better estimate. Now let's use this approach to obtain an envelope of speech signal.

Speech signal

from scipy.io import wavfile # package to read WAV file

# read WAV file and plot data
[fs, x] = wavfile.read('./my_sounds/count.wav')
x = x/np.max(np.abs(x))
t = np.arange(0,len(x))/fs
# get analytic signal
ax = signal.hilbert(x)
envelope_hat = np.abs(ax)
taps = signal.firwin(N, fL_hz/nyq_rate, window=('kaiser', beta))
envelope_hat_filt = signal.filtfilt(taps, 1,envelope_hat)

fig = plt.subplots(figsize=(16,4))
ax = plt.subplot(1,2,1)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope_hat,color='red',label='ENVP. EST.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
# plt.xlim(0,0.5)
plt.ylim(-1,1.5)


ax = plt.subplot(1,2,2)
ax.plot(t,x,color='black',label='SIGNAL')
ax.plot(t,envelope_hat_filt,color='red',label='ENVP. EST. FILT.')
ax.set_xlabel('TIME [in s]')
ax.set_ylabel('A.U.')
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
# plt.xlim(0,0.5)
plt.ylim(-1,1.5)
plt.show()

/Users/neeks/miniconda2/envs/py36/lib/python3.6/site-packages/ipykernel_launcher.py:4: WavFileWarning: Chunk (non-data) not understood, skipping it.
  after removing the cwd from sys.path.

Before we sum-up complete lets see how good is our filter. We will apply it on white noise and see the resulting spectrum before and after application.

# plot signal and spectrum
x = np.random.normal(0,1,size=num_samples)
y = signal.filtfilt(taps, 1,x)
fig  = plt.subplots(figsize=(16,4))

ax = plt.subplot(1,2,1)
ax.plot(np.arange(0,num_samples//2+1,1)/num_samples*fs,10*np.log(abs(np.fft.rfft(x))),label='SIGNAL')
ax.plot(np.arange(0,num_samples//2+1,1)/num_samples*fs,10*np.log(abs(np.fft.rfft(y))),label='SIGNAL FILTERED')
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel('FREQUENCY [in Hz]',fontsize=14)
ax.set_ylabel('SPECTRUM [in dB]',fontsize=14)
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.ylim(10,80)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)

ax = plt.subplot(1,2,2)
ax.plot(taps)
ax.grid(True)
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_xlabel('taps',fontsize=14)
ax.set_ylabel('A.U.',fontsize=14)
ax.legend(loc='upper right',frameon=False,fontsize=12)
plt.xticks(fontsize=13)
plt.yticks(fontsize=13)

plt.show()

No handles with labels found to put in legend.

From the plot it is clear that FIR filter does a good job in removing the spectrum beyond 50 Hz. Lets visualize the filter response and spectrum.

That's it!

To sum it up, assuming you went through the above, we now understand that

Hilbert transform can be used to estimate signal envelope
The estimation is very accurate for tone signals. In general, it is accurate for narrowband carrier signals.
The estimation performance degrades for wideband carriers, like sum of tones or white noise (broadband signal).
The estimation performance improves on applying a lowpass filter to the envelope estimate.
Such lowpass filtering of the envelope estimate can be also applied to speech signals.

The above observations can be reasoned from a theoretical angle. I will try to do it another post.