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

发布时间:2024年01月22日

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

文章来源:https://blog.csdn.net/weixin_39402231/article/details/135739226
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。