Hilbert 变换提取信号特征的 Python 实现 - 文章 - 开发者社区

picture.image

希尔伯特变换(hilbert transform) 一个连续时间信号s(t)的希尔伯特变换等于该信号通过具有冲激响应h(t)=1/πt的线性系统以后的输出响应sh(t)。

好的，这是Hilbert变换的定义，我们这里讨论它的一个具体用途，提取信号特征值，提取信号特征值有什么用呢？

先来一段特征值的定义：

设 A 是n阶方阵，如果存在数m和非零n维列向量 x，使得 Ax=mx 成立，则称 m 是A的一个特征值或本征值。

所以信号特征值可以用来代替一段信号对系统产生的影响，或者说取代一段信号并填补它功能的空缺，我们按照这个思想，在外界或内部改变信号的某些条件之后，可以用特征值随外界或内部因素的变化图像来反映干扰因素对信号的影响，或者反映信号本身对外界干扰的抵抗力及自身的稳定性、鲁棒性等性质。

那么具体怎么用呢？设有一个实值函数s(t)，其希尔伯特反变换记作 \hat{s}(t) ，则有：

希尔伯特变换：

picture.image

希尔伯特反变换：

picture.image

上面说的就是初始信号 s(t) ，我们首先把它做一次希尔伯特变换，然后提取特征值。

都有什么特征值呢？

那我们得先补充一个知识点，包络：随机过程的振幅随着时间变化的曲线。也就是信号的振幅与时间 （A-t） 间的函数关系。

那么根据不同信号包络上高阶统计量特征不同，可分为R值特征和J值特征。

信号包络R值特征：R值是信号包络的方差与包络均值的平方的比值：

picture.image

信号包络J值特征：J值是信号包洛的四阶矩和二阶矩平方之间的差值与包络均值的平方的比值：

picture.image

*双谱特征：利用积分的方法将二维双谱积分函数转换为一维双谱积分函数，从而实现计算方法简单的积分双谱：

picture.image

根据辐射源信号指纹识别技术 http://xueshu.baidu.com/usercenter/paper/show?paperid=f2c7b987bd3af320909da0e1c09a723b&site=xueshu\_se ，这篇论文，三种特征值能够很好的表现信号性能。

好了，基础的知识就是这些，下面我们写代码实现Hilbert变换，计算并提取三种特征值。


          
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)    #产生一个给定形状和类型的用0填充的数组   
  
for i in range(len(t)):   
  
    m[i] = a[int(math.floor(t[i]))]   
  
ts1 = np.arange(0, (np.int64(1/sampling_t) * size))/(10*(m+1))   
  
fsk = np.cos(np.dot(2 * pi, ts1) + pi / 4)   
  
def awgn(y, snr):      #snr为信噪比dB值   
  
    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):        #[feature1, f2, f3] = rj(noise\_bpsk, fs)   
  
    h = fftpack.hilbert(y)  # hilbert变换   
  
    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  # 行数10   
  
    nrecs = np.int64(1 / sampling_t)  # 列数100   
  
    nlag = 20   
  
    nsamp = nrecs  # 每段样本数100   
  
    nrecord = size   
  
    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(nsamp)   
  
    for k in range(nrecord):   
  
        x = y[k][ind]   
  
        x = x - np.mean(x)   
  
        for j in range(nlag + 1):   
  
            z = np.multiply(x[np.arange(nsamp - j)], x[np.arange(j, nsamp)])   
  
            for i in range(j, nlag + 1):   
  
                sum = np.mat(z[np.arange(nsamp - i)]) * np.mat(x[np.arange(i, nsamp)]).T   
  
                sum = sum / nsamp   
  
                c3[i][j] = c3[i][j] + sum  # i,j顺序   
  
    c3 = c3 / nrecord   
  
    c3 = c3 + np.mat(np.tril(c3, -1)).T  # 取对角线以下三角,c3为矩阵   
  
    c31 = c3[1:, 1:]   
  
    c32 = np.mat(np.zeros((nlag, nlag), dtype=np.float32))   
  
    c33 = np.mat(np.zeros((nlag, nlag), dtype=np.float32))  # 不可以直接3者相等   
  
    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:, 0])  # 上下翻转,翻转后依然为矩阵   
  
            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, :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))))          #41*41   
  
    Bspec = fftpack.fft2(cmat, [nfft, nfft])       
  
    Bspec = np.fft.fftshift(Bspec)                #128*128   
  
    waxis = np.arange(-nfft / 2, nfft / 2) / nfft   
  
    X, Y = np.meshgrid(waxis, waxis)   
  
    plt.contourf(X, Y, abs(Bspec),alpha=0,cmap=plt.cm.hot)   
  
    plt.contour(X, Y, abs(Bspec))   
  
    plt.show()   
  
    return Bspec   
  
def features(s):   
  
#    for mc in range(2):   
  
        snr = np.random.uniform(0, 20)      #从一个均匀分布集合中随机采样，左闭右开--[low,high)   
  
        s = awgn(s,snr)            #在原始信号的基础上增加SNR信噪比的噪音   
  
        rj = np.array(feature_rj(s))              #计算R,J特征   
  
        z = feature_Bispectrum(s)                  #计算双谱特征，并画图像   
  
        xx = np.int64(np.sqrt(np.size(z))/2)   
  
        z = np.array(z[:xx,xx:])   
  
        z = np.tril(z).real              #取复数z的实部   
  
        bis = np.zeros((1, xx),dtype=np.float32)    #零组   
  
        for i in range(xx):   
  
            for j in range(xx-i):   
  
                bis[0][i] = bis[0][i]+z[xx-1-j][i+j]   
  
        m = bis[0].reshape(1,xx)   
  
        normalized = preprocessing.normalize(m)[0,:]    #样本各个特征值除以各个特征值的平方之和   
  
        features = np.hstack((rj,normalized))          #合并数组r,j和normalized   
  
        return features   
  
print(features(fsk))

好的，特征值已经有了，那么下面我们怎么使用这些特征值来描述初始信号的变化趋势呢？

画出特征值曲线是个不错的办法，趋势一目了然，为了说明信号受到影响，我们分别从信号的振幅 A ，角频率 \omega ，初始频率 S 的改变来说明问题。初始信号图像如下：

picture.image

为了统一起见，我们分别令 A ， \omega ， S 从0到 2\Pi 均匀变化，对于振幅 A 的变化，补充代码计算 R ， J 特征值，画出双谱图：


          
R = []  
  
J = []   
  
ts1 = np.arange(0, (np.int64(1/sampling_t) * size))/(10*(m+1))   
  
W = []   
  
Z = []   
  
for i in range(0,40,1):   
  
    W.append(i / (2 * pi))   
  
for a1 in W:   
  
    global r,j   
  
    fsk = a1 * np.cos(np.dot(2 * pi, ts1) + pi / 4)   
  
    features(fsk)   
  
    R.append(r)   
  
    J.append(j)   
  
plt.plot(W, R, color='green', label='1')   
  
plt.legend() # 显示图例   
  
plt.xlabel('A[0-2 * pi]')   
  
plt.ylabel('R')   
  
plt.show()   
  
plt.plot(W, J, color='red', label='2')   
  
plt.legend() # 显示图例   
  
plt.xlabel('A[0-2 * pi]')   
  
plt.ylabel('J')   
  
plt.show()   
  
plt.plot(W, Z, color='red', label='3')   
  
plt.legend() # 显示图例   
  
plt.xlabel('A[0-2 * pi]')   
  
plt.ylabel('trend')   
  
plt.show()

R 特征值变化图像：

picture.image

J 特征值变化图像：

picture.image

双谱图像：

picture.image

对于角频率 ω 的变化，补充代码计算 R ， J 特征值，画出 R 、 J 、双谱图像：


          
R = []  
  
J = []   
  
ts1 = np.arange(0, (np.int64(1/sampling_t) * size))/(10*(m+1))   
  
W = []   
  
Z = []   
  
for i in range(0,40,1):   
  
    W.append(i / (2 * pi))   
  
for w in W:   
  
    global r,j   
  
    fsk = np.cos(np.dot(w, ts1) + pi / 4)   
  
    features(fsk)   
  
    R.append(r)   
  
    J.append(j)   
  
plt.plot(W, R, color='green', label='1')   
  
plt.legend() # 显示图例   
  
plt.xlabel('w[0-2 * pi]')   
  
plt.ylabel('R')   
  
plt.show()   
  
plt.plot(W, J, color='red', label='2')   
  
plt.legend() # 显示图例   
  
plt.xlabel('w[0-2 * pi]')   
  
plt.ylabel('J')   
  
plt.show()   
  
plt.plot(W, Z, color='red', label='3')   
  
plt.legend() # 显示图例   
  
plt.xlabel('A[0-2 * pi]')   
  
plt.ylabel('trend')   
  
plt.show()

R 特征值变化图像：

picture.image

J 特征值变化图像：

picture.image

双谱图像：

picture.image

对于初始频率 S 的变化，补充代码计算 R ， J 特征值，画出双谱图：


          
R = []  
  
J = []   
  
ts1 = np.arange(0, (np.int64(1/sampling_t) * size))/(10*(m+1))   
  
W = []   
  
Z = []   
  
for i in range(0,40,1):   
  
    W.append(i / (2 * pi))   
  
for s in W:   
  
    global r,j   
  
    fsk = np.cos(np.dot(2 * pi, ts1) + s)   
  
    features(fsk)   
  
    R.append(r)   
  
    J.append(j)   
  
plt.plot(W, R, color='green', label='1')   
  
plt.legend() # 显示图例   
  
plt.xlabel('s[0-2 * pi]')   
  
plt.ylabel('R')   
  
plt.show()   
  
plt.plot(W, J, color='red', label='2')   
  
plt.legend() # 显示图例   
  
plt.xlabel('s[0-2 * pi]')   
  
plt.ylabel('J')   
  
plt.show()   
  
plt.plot(W, Z, color='red', label='3')   
  
plt.legend() # 显示图例   
  
plt.xlabel('A[0-2 * pi]')   
  
plt.ylabel('trend')   
  
plt.show()

R 特征值变化图像：

picture.image

J 特征值变化图像：

picture.image

双谱图像：

picture.image

二维的双谱图像看起来不是很方便，那么对图像进行降维处理，像这种比较对称的矩阵，用矩阵平均值是最适合的了，这里只需将希尔伯特变换后的包络矩阵求下均值，改变下述代码即可：


          
Bspec = fftpack.fft2(cmat, [nfft, nfft])  
  
Bspec = np.fft.fftshift(Bspec)                #128*128   
  
waxis = np.arange(-nfft / 2, nfft / 2) / nfft   
  
X, Y = np.meshgrid(waxis, waxis)   
  
plt.contourf(X, Y, abs(Bspec))   
  
plt.contour(X, Y, abs(Bspec))   
  
Z.append(np.mean(abs(Bspec)))

则双谱特征降维后图像如下，振幅A的双谱特征变化情况：

picture.image