Python环境下一维时间序列信号的时频脊线追踪方法

瞬时频率是分析调频信号的一个重要参数，它表示了信号中的特征频率随时间的变化。使用短时傅里叶变换或小波变换获得信号的时频表示TFR后，从TFR中估计信号各分量的瞬时频率，即可获得信号中的特征信息。在TFR中，调频信号的特征分量通常显示为一条条线性或非线性变化的能量带，虽然能量带在TFR中清晰可见，然而无法通过肉眼观察从中估计信号特征分量的瞬时频率。在能量带中每个时刻上能量最高的点，被称为时频脊线。时频脊线反映了特征分量的频率随时间的变化情况，因此通过提取TFR中的时频脊线可估计调频信号特征分量的瞬时频率。

下面给出几个小demo，首选加载相关模块

import numpy as np
import pywt
from scipy.signal import chirp,sweep_poly
import matplotlib.pyplot as plt
from tfridge import ridge as rt

定义可视化函数

def visualize(signal, tf_transf,ridge,flip_plot=False):
    plt.plot(signal)
    plt.show()
    
    plt.figure()
    if(flip_plot):
        plt.imshow(np.flipud(np.abs(tf_transf)), aspect='auto', vmin=0, vmax=np.max(np.abs(tf_transf)), cmap='jet')
        plt.plot(len(tf_transf)-ridge,linestyle = '--',color='black')
    else:
        plt.imshow(np.abs(tf_transf), aspect='auto', cmap='jet')
        plt.plot(ridge,linestyle = '--',color='black')
    
    plt.show()

第一个简单的例子：

test_matrix=np.array([[1, 4, 4],[2, 2, 2],[5,5,4]])
fs_test =np.exp(np.array([1,2,3]))

Energy,ridge_idx,_ = rt.extract_fridges(test_matrix,fs_test,penalty=2.0)

print('ridge follows indexes:', ridge_idx)

ridge follows indexes: [[2]

[2]

[2]]

第2个简单的例子

###  Two constant frequency signals
    
sig_lgth = 500
t_vec=np.linspace(0,10,sig_lgth)
f1=0.5
f2=2.0
signal_1=np.sin(f1*2*np.pi*t_vec)
signal_2=np.cos(f2*2*np.pi*t_vec)

signal=signal_1+signal_2

nbrVoices=32
nbr_scales=279
scales=np.exp(np.linspace(np.log(1.51),np.log(622.207),nbr_scales))
cwtmatr, freqs = pywt.cwt(signal,scales,'cmor2.0-1.0')

penalty=2.0
# CWT example
Energy,ridge_idx,_ = rt.extract_fridges(cwtmatr,scales,penalty,num_ridges=2,BW=25)
plt.figure()
visualize(signal, cwtmatr, ridge_idx)

第3个简单的例子：

### Two chirp signals with linear and quadratic frequrncy variation

sign_chirp_1 = chirp(t_vec, f0=2, f1=10, t1=20, method='linear')
sign_chirp_2 = chirp(t_vec, f0=.4, f1=7, t1=20, method='quadratic')

sign_chirp=sign_chirp_1+sign_chirp_2


scales=np.exp(np.linspace(np.log(1.51),np.log(622.207),nbr_scales))
cwtmatr, freqs = pywt.cwt(sign_chirp,scales,'cmor2.0-1.0')

penalty=.3

# CWT example
Energy,ridge_idx,_ = rt.extract_fridges(cwtmatr,scales,penalty,num_ridges=2,BW=25)
plt.figure()
visualize(sign_chirp, cwtmatr, ridge_idx)

第4个简单的例子：

####  One sweep signal and one constant frequency signal


p1 = np.poly1d([0.025, -0.36, 1.25, 2.0])

sweep_sig = sweep_poly(t_vec, p1)+signal_1

scales=np.exp(np.linspace(np.log(1.51),np.log(622.207),nbr_scales))
cwtmatr, freqs = pywt.cwt(sweep_sig,scales,'cmor2.0-1.5')

penalty=2.0

# CWT example
Energy,ridge_idx,_ = rt.extract_fridges(cwtmatr,scales,penalty,num_ridges=2,BW=25)
plt.figure()
visualize(sweep_sig, cwtmatr, ridge_idx)

工学博士，担任《Mechanical System and Signal Processing》审稿专家，担任
《中国电机工程学报》优秀审稿专家，《控制与决策》，《系统工程与电子技术》，《电力系统保护与控制》，《宇航学报》等EI期刊审稿专家，担任《计算机科学》，《电子器件》 ， 《现代制造过程》 ，《电源学报》，《船舶工程》 ，《轴承》 ，《工矿自动化》 ，《重庆理工大学学报》 ，《噪声与振动控制》 ，《机械传动》 ，《机械强度》 ，《机械科学与技术》 ，《机床与液压》，《声学技术》，《应用声学》，《石油机械》，《西安工业大学学报》等中文核心审稿专家。
擅长领域：现代信号处理，机器学习，深度学习，数字孪生，时间序列分析，设备缺陷检测、设备异常检测、设备智能故障诊断与健康管理PHM等。