Python的matplotlib绘制小波能量谱

程序员文章站 2024-03-26 12:41:53

...

时间小波能量谱

反映信号的小波能量沿时间轴的分布。

由于小波变换具有等距效应，所以有：

∫ R ∣ f ( t ) ∣ 2 d x = C ψ − 1 ∫ R ∫ R ∣ W T f ( a , b ) ∣ 2 d a d b a 2 \int_{R}|f(t)|^{2}{\rm d}x=C^{-1}_{\psi}\int_{R}\int_{R}|WT_{f}(a,b)|^{2}\frac{{\rm d}a{\rm d}b}{a^2} ∫R∣f(t)∣2dx=Cψ−1∫R∫R∣WTf(a,b)∣2a2dadb ①

式中 ∣ W T f ( a , b ) ∣ 2 |WT_{f}(a,b)|^2 ∣WTf(a,b)∣2 表示信号尺度，对式①在平移因子b方向上进行加权积分：

∫ − ∞ + ∞ ∣ f ( t ) ∣ 2 d x = ∫ − ∞ + ∞ E ( b ) d b \int^{+\infty}_{-\infty}|f(t)|^{2}{\rm d}x=\int^{+\infty}_{-\infty}E(b){\rm d}b ∫−∞+∞∣f(t)∣2dx=∫−∞+∞E(b)db ②

式中： E ( b ) = C ψ − 1 ∫ − ∞ + ∞ ∣ W T f ( a , b ) ∣ 2 a 2 d a E(b)=C^{-1}_{\psi}\int^{+\infty}_{-\infty}|WT_{f}(a,b)|^2\frac{a^{2}}{{\rm d}a} E(b)=Cψ−1∫−∞+∞∣WTf(a,b)∣2daa2 代表时间-小波能量谱。

尺度小波能量谱

反映信号的小波能量随尺度的变化情况。

同理，对式①在尺度方向上进行加权积分：

∫ − ∞ + ∞ ∣ f ( t ) ∣ 2 d t = C ψ − 1 ∫ − ∞ + ∞ a − 2 E ( a ) d a \int^{+\infty}_{-\infty}|f(t)|^{2}{\rm d}t=C^{-1}_{\psi}\int^{+\infty}_{-\infty}a^{-2}E(a){\rm d}a ∫−∞+∞∣f(t)∣2dt=Cψ−1∫−∞+∞a−2E(a)da ③

式中： E ( a ) = ∫ − ∞ + ∞ ∣ W T f ( a , b ) ∣ 2 d b E(a)=\int^{+\infty}_{-\infty}|WT_{f}(a,b)|^2{\rm d}b E(a)=∫−∞+∞∣WTf(a,b)∣2db 代表尺度小波能量谱。

连续小波变换

连续小波变换的结果是一个小波系数矩阵，随着尺度因子和位移因子变化。然后将系数平方后得到小波能量，把每个尺度因子对应的所有小波能量进行叠加，那么就可以得到随尺度因子变换的小波能量谱曲线。把尺度换算成频率后，这条曲线就可视为是频谱图。

代码如下：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import pywt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import MultipleLocator, FormatStrFormatter

# 解决负号显示问题
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
plt.rcParams.update({'text.usetex': False, 'font.family': 'serif', 'font.serif': 'cmr10', 'mathtext.fontset': 'cm'})

font1 = {'family': 'SimHei', 'weight': 'normal', 'size': 12}
font2 = {'family': 'Times New Roman', 'weight': 'normal', 'size': 18}
label = {'family': 'SimHei', 'weight': 'normal', 'size': 15}

xlsx_path = "../小波能量谱/小波能量谱作图.xlsx"
sheet_name = "表名"      
data_arr = pd.read_excel(xlsx_path, sheet_name=sheet_name)

column_name = '信号数据'    
row = 1024
y = data_arr[column_name][0:row]
x = data_arr['里程'][0:row]

scale = np.arange(1, 50)
wavelet = 'gaus1'   # 'morl'  'gaus1'  小波基函数

# 时间-尺度小波能量谱
def time_scale_spectrum():
    coefs, freqs = pywt.cwt(y, scale, wavelet)  # np.arange(1, 31) 第一个参数必须 >=1     'morl'  'gaus1'
    scale_freqs = np.power(freqs, -1)  # 对频率freqs 取倒数变为尺度

    fig = plt.figure(figsize=(5, 4))
    ax = Axes3D(fig)

    # X:time   Y:Scale   Z:Amplitude
    X = np.arange(0, row, 1)  # [0-1023]
    Y = scale_freqs
    X, Y = np.meshgrid(X, Y)
    Z = abs(coefs)

    # 绘制三维曲面图
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')

    # 设置三个坐标轴信息
    ax.set_xlabel('$Mileage/km$', color='b', fontsize=12)
    ax.set_ylabel('$Scale$', color='g', fontsize=12)
    ax.set_zlabel('$Amplitude/mm$', color='r', fontsize=12)

    plt.draw()
    plt.show()


# 时间小波能量谱
def time_spectrum():
    coefs, freqs = pywt.cwt(y, scale, wavelet)
    coefs_pow = np.power(coefs, 2)      # 对二维数组中的数平方
    spectrum_value = [0] * row    # len(freqs)

    # 将二维数组按照里程叠加每个里程上的所有scale值
    for i in range(row):
        sum = 0
        for j in range(len(freqs)):
            sum += coefs_pow[j][i]
        spectrum_value[i] = sum

    fig = plt.figure(figsize=(7, 2))

    line_width = 1
    line_color = 'dodgerblue'
    line_style = '-'

    T1 = fig.add_subplot(1, 1, 1)
    T1.plot(x, spectrum_value, label='模拟', linewidth=line_width, color=line_color, linestyle=line_style)
    T1.legend(loc='upper right', prop=font1, frameon=True)  # lower ,left

    # 坐标轴名称
    T1.set_xlabel('里程/$km$', fontsize=15, fontdict=font1)  # fontdict设置子图字体
    T1.set_ylabel('$E/mm^2$', fontsize=15, fontdict=font1)
    # 坐标刻度值字体大小
    T1.tick_params(labelsize=15)

    print(spectrum_value[269])
    plt.show()


# 尺度小波能量谱
def scale_spectrum():
    coefs, freqs = pywt.cwt(y, scale, wavelet)
    coefs_pow = np.power(coefs, 2)      # 对二维数组中的数平方
    scale_freqs = np.power(freqs, -1)   # 对频率freqs 取倒数变为尺度
    spectrum_value = [0] * len(freqs)    # len(freqs)

    # 将二维数组按照里程叠加每个里程上的所有scale值
    for i in range(len(freqs)):
        sum = 0
        for j in range(row):
            sum += coefs_pow[i][j]
        spectrum_value[i] = sum

    fig = plt.figure(figsize=(7, 4))

    line_width = 1
    line_color1 = 'dodgerblue'
    line_style1 = '-'

    T1 = fig.add_subplot(1, 1, 1)
    T1.plot(scale_freqs, spectrum_value, label=column_name, linewidth=line_width, color=line_color1, linestyle=line_style1)
    T1.legend(loc='upper right', prop=font1, frameon=True)  # lower ,left

    # 坐标轴名称
    T1.set_xlabel('$Scale/mm^2$', fontsize=15, fontdict=font1)  # fontdict设置子图字体
    T1.set_ylabel('$E/mm^2$', fontsize=15, fontdict=font1)
    # 坐标刻度值字体大小
    T1.tick_params(labelsize=15)

    plt.show()


# 通过调用下面三个不同的函数选择绘制能量谱
# time_scale_spectrum()
# time_spectrum()
# scale_spectrum()