连续小波变换python代码实现:
import matplotlib.pyplot as plt
import numpy as np
import pywt
def MyCWT(y, fs, wavelet='cmor2-5', total_scal=512):
'''
连续小波变换CWT
:param y: 信号,nnumpy, (n,)
:param fs: 采样频率
:param wavelet: 复小波cmorB-C,推荐‘cmor3-3’,‘cmor3-5’
:param total_scal: 尺度
:return:
t -- 时间轴
f -- 频率轴
z -- 幅值轴
Notes:
total_scal越大,频率分辨率越高,但是计算消耗越大,一般推荐256/512
'''
N = len(y)
t = np.linspace(0, (N-1)/fs, N) # 时间轴
fc = pywt.central_frequency(wavelet) # 小波中心频率
cparam = 2. * fc * total_scal
scales = cparam / np.arange(1, total_scal + 1) # 各个尺度,用于获得频率序列
coefficients, f = pywt.cwt(y, scales=scales, wavelet=wavelet, sampling_period=1.0 / fs)
z = np.abs(coefficients)
return t, f, z
if __name__ == "__main__":
fs = 12800 # 采样频率
f0 = 60 # 信号基频
f1 = 90 # 信号基频
t_total = 0.32 # 信号持续时间
t = np.linspace(0, t_total, int(fs * t_total), endpoint=False)
signal = np.sin(2 * np.pi * f0 * t) + np.sin(2 * np.pi * f1 * t)
t, f, z = MyCWT(signal, fs)
plt.figure(figsize=(6, 4), constrained_layout=True)
plt.contourf(t, f, z, cmap='jet')
plt.xlabel("Time (s)")
plt.ylabel("Frequency (Hz)")
plt.ylim(0, 500) # 限制频率范围
plt.colorbar(label="Magnitude")
plt.show()
结果如下: