Extracting Signal Features with Hilbert Transform and Bispectrum in Python

The Hilbert transform of a continuous‑time signal s(t) is the output of a linear system with impulse response h(t)=1/πt, denoted sh(t). Using this transform we can extract signal feature values that replace the signal’s effect on a system and reveal interference, stability, and robustness.

For a real‑valued function s(t) we denote its Hilbert inverse as \hat{s}(t). The Hilbert transform and its inverse are illustrated below:

Signal envelope (A‑t) is the amplitude curve of a random process. Two envelope‑based statistics are defined:

R value : the variance of the envelope divided by the square of its mean.

J value : the difference between the fourth‑order moment and the square of the second‑order moment, divided by the square of the mean.

These are shown in the following formulas (images):

The bispectrum feature is obtained by converting the two‑dimensional bispectrum integral into a one‑dimensional form, simplifying computation.

import numpy as np
from math import pi
import matplotlib.pyplot as plt
import math
from scipy import fftpack
from sklearn import preprocessing
import neurolab as nl

size = 10
sampling_t = 0.01
t = np.arange(0, size, sampling_t)
a = np.random.randint(0, 2, size)
m = np.zeros(len(t), dtype=np.float32)
for i in range(len(t)):
    m[i] = a[int(math.floor(t[i]))]

ts1 = np.arange(0, (int(1/sampling_t) * size)))/(10*(m+1))
fsk = np.cos(2*pi*ts1 + pi/4)

def awgn(y, snr):
    snr = 10**(snr/10.0)
    xpower = np.sum(y**2)/len(y)
    npower = xpower / snr
    return np.random.randn(len(y))*np.sqrt(npower) + y

def feature_rj(y):
    h = fftpack.hilbert(y)
    z = np.sqrt(y**2 + h**2)
    m2 = np.mean(z**2)
    m4 = np.mean(z**4)
    r = abs((m4-m2**2)/m2**2)
    Ps = np.mean(y**2)/2
    j = abs((m4-2*m2**2)/(4*Ps**2))
    return (r, j)

def feature_Bispectrum(y):
    ly = size
    nrecs = int(1/sampling_t)
    nlag = 20
    nfft = 128
    Bspec = np.zeros((nfft, nfft), dtype=np.float32)
    y = y.reshape(ly, nrecs)
    c3 = np.zeros((nlag+1, nlag+1), dtype=np.float32)
    ind = np.arange(nrecs)
    for k in range(ly):
        x = y[k][ind] - np.mean(y[k][ind])
        for j in range(nlag+1):
            z = np.multiply(x[np.arange(nrecs-j)], x[np.arange(j, nrecs)])
            for i in range(j, nlag+1):
                sum = np.mat(z[np.arange(nrecs-i)]) * np.mat(x[np.arange(i, nrecs)]).T
                sum = sum / nrecs
                c3[i][j] += sum
    c3 = c3 / ly
    c3 = c3 + np.mat(np.tril(c3, -1)).T
    c31 = c3[1:,1:]
    c32 = np.mat(np.zeros((nlag, nlag), dtype=np.float32))
    c33 = np.mat(np.zeros((nlag, nlag), dtype=np.float32))
    c34 = np.mat(np.zeros((nlag, nlag), dtype=np.float32))
    for i in range(nlag):
        x = c31[i:, i]
        c32[nlag-1-i, 0:nlag-i] = x.T
        c34[0:nlag-i, nlag-1-i] = x
        if i < (nlag-1):
            x = np.flipud(x[1:])
            c33 = c33 + np.diag(np.array(x)[:,0], i+1) + np.diag(np.array(x)[:,0], -(i+1))
    c33 = c33 + np.diag(np.array(c3)[0,:-1])
    cmat = np.vstack((np.hstack((c33, c32, np.zeros((nlag,1), dtype=np.float32))),
                     np.hstack((np.vstack((c34, np.zeros((1,nlag), dtype=np.float32))), c3))))
    Bspec = fftpack.fft2(cmat, [nfft, nfft])
    Bspec = np.fft.fftshift(Bspec)
    return Bspec

def features(s):
    snr = np.random.uniform(0, 20)
    s = awgn(s, snr)
    r, j = feature_rj(s)
    B = feature_Bispectrum(s)
    return np.hstack((r, j, np.mean(np.abs(B))))

print(features(fsk))

The script then varies amplitude A, angular frequency ω, and initial frequency S from 0 to 2π, computes R, J, and bispectrum features for each case, and plots the resulting curves. Sample plots for R, J, and bispectrum versus A, ω, and S are shown below:

Because the two‑dimensional bispectrum matrix is symmetric, its dimensionality can be reduced by averaging the matrix values. The modified code replaces the bispectrum visualization with a single scalar:

Bspec = fftpack.fft2(cmat, [nfft, nfft])
Bspec = np.fft.fftshift(Bspec)
Z.append(np.mean(abs(Bspec)))

Resulting reduced‑dimensional bispectrum plots for varying amplitude, angular frequency, and initial frequency are displayed respectively: