瞬时频率是分析调频信号的一个重要参数,它表示了信号中的特征频率随时间的变化。使用短时傅里叶变换或小波变换获得信号的时频表示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等。