当前位置：首页 > news >正文

python代码实现离散haar小波变换和db4小波变换

news 2025/4/2 21:09:13

1.小波变换

傅里叶变换

短时傅里叶变换（STFT）

连续小波变换（CWT）

离散小波变换（DWT）

重构信号

2.快速算法实现（Mallat算法）

haar小波

1.分解代码

2.重构代码

3.优化代码，保留多级分解系数

4.音频信号测试

Daubechies(dbN)小波系

1.介绍

2.分解代码

3.重构代码

4.优化代码

3.真实信号测试

1.信号读取

2.信号处理（分解，滤波，重构）

3.信号分析对比（频谱分析，信噪比等）

1.小波变换

小波变换是指，输入信号通过和小波函数做卷积，被分解为一系列的信号，对应着不同的时域信息和频域信息。因此，小波变换可以实现时频分析。

对信号的分解和重构需要用到一组正交基，在傅里叶变换当中是不同频率的正弦波，而在小波变换中是一系列的小波函数。

傅里叶变换

傅里叶变换用不同频率的正弦波作为基，变换过程如下：

傅里叶变换的缺点是只能显示频率成分而不能确切知道这些频率成分出现的时间，有三个不同的信号，但他们的频谱图相差无几，这样仅凭频谱图无法将信号还原。

平稳信号是指分布参数或者分布律随时间不发生变化的信号。也就是说，平稳信号的统计特性不随时间变化而变化。

非平稳信号是指分布参数或者分布律随时间发生变化的信号。也就是说，非平稳随机信号的统计特征是时间的函数（随时间变化），如下图所示：

平稳信号	非平稳信号

对此采用加窗傅里叶变换，来处理非平稳信号。

短时傅里叶变换（STFT）

为了弥补FT的不足，把整个时间域分解成无数个等长的小过程（加窗），每个过程近似平稳，再傅里叶变换，就知道在哪个时间点上出现了什么频率。

新的问题是对于窗宽度的选取，窗太窄，窗内的信号太短，会导致频率分析不够精准，频率分辨率差。窗太宽，时域上又不够精细，时间分辨率低。

对于时变的非稳态信号，高频适合小窗口，低频适合大窗口。然而STFT的窗口是固定的，在一次STFT中宽度不会变化，所以STFT还是无法满足非稳态信号变化的频率的需求。

因此采用小波变换，小波变换与傅里叶变换的区别就是选取的基不同，傅里叶变换分解得到了一系列不同频率的正弦波，小波变换得到了一系列不同尺度的小波。

连续小波变换（CWT）

把小波函数想象成一个观察者，信号与之卷积，得到的结果则代表了信号中含有该成分的多少，这一系列的观察者可以称为一组基，用一组正交基就能够实现信号的分解与重构。

可以参考b站视频。

【一次讲透小波变换原理，全新角度切入，20分钟时长警告！】 https://www.bilibili.com/video/BV1tu4y1M7th/?share_source=copy_web&vd_source=2988ac044a7ec04a360fa2edbb97cefe

其中还讲到了如何构造滤波器系数等。

连续小波变换的公式如下：

取不同的尺度参数可以对小波函数进行尺度变换，即伸缩。而不同的平移参数可以确定不同的时间，这样就能做到时频分析，确定不同的频率在什么时间出现。

用不同尺度的小波函数分别对信号进行卷积

把小波w(t)和原函数f(t)的开始部分进行比较（实际上就是作內积），计算系数C。系数C表示该部分函数与小波的相似程度。
把小波向右移k单位，得到小波w(t-k)，重复1。重复该步骤直至函数f结束.
扩展小波w(t)，得到小波w(t/2)，重复步骤1,2.
不断扩展小波，重复1,2,3

变换后可以得到如下的标度图，两个坐标轴分别为time和Scales（即b和a）Lower Scales对应高频成分，High Scales对应低频成分。

一些知名的小波基函数

下面是一个基于时间对比，频率对比，以及短时间傅立叶分析的对比：

海森堡不确定性

由于输入的序列信号是离散的，因此讨论离散情况下如何进行小波变换。

离散小波变换（DWT）

DWT 通过离散化的尺度（scale）和平移（translation）参数，将信号分解为不同频带的近似系数（Approximation Coefficients）和细节系数（Detail Coefficients）。其核心公式如下：

小波变换有母小波和父小波：

母小波是小波函数 $\psi$ （t），做高通滤波保留细节；

父小波为尺度函数 $\varphi$ （t），做低通滤波保留近似形状；

如图经过离散变化将会得到两个分量，原信号为带噪声的正弦波序列，

而分离出来的approximation代表低通滤波得到的信号包络，也就是大致的轮廓

detail则反映了高通滤波得到的细节

以haar小波为例

若进行多级分解，则可将上一级的近似部分再分解为近似部分和细节部分。

和连续小波变换一样，也会得到标度图，下图是DWT的标度图

重构信号

通过近似系数和细节系数的线性组合，可完全重建原始信号 x[n]x[n]。

噪声一般会在高频部分，对高频部分进行滤波后再重构信号，会有去噪效果。

如图，其中CD1和CD2为高频部分，将其置零后再重构信号。

一般重构的方法会用到阈值限定，在一定阈值下对高频部分进行滤波，因为高频部分也会有我们需要的信息，阈值限定大致分为硬阈值和软阈值。

硬阈值（Hard Thresholding）：定下一个阈值T，大于T的点全部置零，缺点是会引入信号不连续现象。

软阈值（Soft Thresholding）：对大于阈值的系数进行收缩（向零压缩），结果更平滑，但可能过度衰减信号幅度。

阈值的选择：

介绍一下通用阈值和分层阈值，后面可能会用到

通用阈值：

分层阈值：

注意：为了计算更方便，本文标准差不使用中位数而改用信号均值代替

2.快速算法实现（Mallat算法）

实际计算中，DWT 通过多分辨率分析（MRA）和滤波器组高效实现：

在多分辨分析中，如正交小波变换可以等效为一组镜像滤波的过程，即信号通过一个分解高通滤波器和分解低通滤波器，自然的高通滤波器输出对应的信号的高频分量部分，称为细节分量，低通滤波器输出对应了信号的相对较低的频率分量部分，称为近似分量。

对应的快速算法称为Mallat算法。

滤波器系数可以通过以下代码获取

import pywt#获取小波变换db4小波的滤波器系数
if __name__ == '__main__':h = pywt.Wavelet('db4').dec_lo  # 低通g = pywt.Wavelet('db4').dec_hi  # 高通print("低通 h:", h)print("高通 g:", g)

得到的结果为

低通 h: [-0.010597401785069032, 0.0328830116668852, 0.030841381835560764, -0.18703481171909309, -0.027983769416859854, 0.6308807679298589, 0.7148465705529157, 0.2303778133088965]
高通 g: [-0.2303778133088965, 0.7148465705529157, -0.6308807679298589, -0.027983769416859854, 0.18703481171909309, 0.030841381835560764, -0.0328830116668852, -0.010597401785069032]

设置好滤波器后直接与序列信号进行卷积即可实现小波变换，遇到的问题是关于序列边界的处理，例如序列长度为n，db4的滤波器长度为8，二者卷积需要给序列补零，

通过高通滤波器和低通滤波器将其分解为高频部分和低频部分，若进行多级分解，则对低频部分重复操作就好。

以haar小波为例，编写DWT。

haar小波

Haar小波是最常用的小波基，公式和图像如下：

实现代码如下

关于代码有一些需要说明的，可以用一些数学库使得代码更简便，我的代码是手搓的，没有使用库，所以会更长一些，如果有需要可以直接让deepseek等给你生成一下就好了。

1.分解代码

import math
import numpy as np
import matplotlib.pyplot as pltdef haarfj(input_signal, t=1, decomposition_level=1):"""对输入信号进行haar小波的多级分解input_signal；输入信号t:信号长度decomposition_level:分解级数return:approx, detail 近似系数和细节系数"""current_signal = input_signal# 构建主循环for level in range(decomposition_level):n = len(current_signal)print(f"\n=== 正在进行第 {level + 1} 级分解 ===")print(f"当前信号长度: {n}")# 判断是否为偶数序列if n % 2 == 0:n = nelse:if n < 2:print("函数无法分解，程序终止")return False  # 返回 False 并跳出函数else:n = n - 1approx = np.zeros(n // 2)  # 当前层近似系数detail = np.zeros(n // 2)  # 当前层细节系数# 计算系数haar小波系数for i in range(n // 2):approx[i] = (current_signal[2 * i] + current_signal[2 * i + 1]) / math.sqrt(2)detail[i] = (current_signal[2 * i] - current_signal[2 * i + 1]) / math.sqrt(2)# 准备下一层分解（使用当前层的近似系数）current_signal = approx# 显示分解结果plt.figure(figsize=(10, 6))plt.subplot(3, 1, 1)plt.title("approx")plt.plot(np.linspace(0, t, len(approx)), approx, label=r'$mxb$')plt.subplot(3, 1, 2)plt.title("detail")plt.plot(np.linspace(0, t, len(detail)), detail, label=r'$mxb$', color='yellow')plt.subplot(3, 1, 3)plt.title("function")plt.plot(np.linspace(0, t, len(input_signal)), input_signal, label=r'$fcb$', color='green')plt.tight_layout()plt.show()return approx, detailif __name__ == '__main__':# 生成示例信号（例如正弦波 + 噪声）t = np.linspace(0, 1, 1000, endpoint=False)noise = 0.1 * np.sin(2 * np.pi * 200 * t)  # 噪声# signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.random.randn(len(t))signal = np.sin(2 * np.pi * 5 * t) + noise#haar小波分解approx, detail = haarfj(signal, 1, 2)

注意：输入的t是信号的长度，不要输入成序列t了，错误示范如下：

如果输入错误，就会得到1000条绘制出来的函数的叠加，下面展示正确的第一级分解：

这是进行第一级分解得到的结果，用这一级的近似系数继续分解，得到第二级，最下方的function展示的一直是原函数，第二级分解如下：

第三级分解：

可以看到，三级分解之后，已经能够滤除大部分噪声。接下来涉及到信号的滤波以及重构部分，先编写重构的代码，要求是输入近似系数和细节系数，同时设置阈值T对细节系数进行滤波，滤波后能够将二者合成上一级的近似系数（或原函数）。采取通用阈值对其进行硬阈值处理，代码如下：

2.重构代码

def haarcg(approx, detail,t = 1):"""进行haar小波变换的重构，重构要对高频成分进行阈值限定从而滤除噪声，本代码只进行一次重构，若有多次，则多次调用即可approx:低频成分detail:高频成分return:reconstructed"""# ============== 1. 参数校验 ==============n = len(detail)if len(approx) != n:raise ValueError("近似系数和细节系数长度必须相同")if n == 0:return []# ============== 2. 计算通用阈值 ==============# 2.1 计算细节系数的均值mean = 0for i in range(n):mean += detail[i] /n# 2.2 计算方差variance = 0for i in range(n):variance += (detail[i] - mean) ** 2 / n# 2.3 计算通用阈值 (Donoho-Johnstone公式)threshold = math.sqrt(variance) * math.sqrt(2 * math.log(n))# ============== 3. 高频系数阈值处理 ==============processed_detail = np.zeros(len(detail))for i in range(len(detail)):# 硬阈值处理：绝对值小于阈值的系数置零if abs(detail[i]) <= threshold:processed_detail[i] = 0else:processed_detail[i] = detail[i]detail = processed_detail# ============== 4. Haar小波重构 ==============reconstructed = np.zeros(2 * n)print(f"=== 正在进行信号重构 ===")for i in range(n):reconstructed[2 * i] = (approx[i] + detail[i]) / math.sqrt(2)if 2 * i + 1 < 2*n:reconstructed[2 * i+1] = (approx[i] - detail[i]) / math.sqrt(2)# ============== 5. 重构结果可视化 ==============plt.figure(figsize=(10, 6))plt.subplot(3, 1, 1)plt.title("approx")plt.plot(np.linspace(0, t, len(approx)), approx, label=r'$mxb$')plt.subplot(3, 1, 2)plt.title("detail")plt.plot(np.linspace(0, t, len(detail)), detail, label=r'$mxb$', color='yellow')plt.subplot(3, 1, 3)plt.title("function")plt.plot(np.linspace(0, t, len(reconstructed)), reconstructed, label=r'$fcb$', color='green')plt.tight_layout()plt.show()return reconstructed

用之前三级分解的系数输入，得到重构的第二级近似系数

重构的第一级

重构的原函数

让噪声更大一些，滤波效果会看的更明显

原信号一级分解	三级分解滤波后重构

左边是第一级分解的图像，右边是三级分解之后又逐级滤波重构后的图像。还可以修改不同的滤波方式，会有不一样的效果。

3.优化代码，保留多级分解系数

由于重构需要用到最后一级的近似系数，以及每一级的细节系数，因此希望在分解时顺便把每一级的细节系数保留下来。

注意1：同时还有需要优化的部分，就是关于序列的长度，分解要求序列长度为偶数，一直分解下去总有奇数的情况，补零或者去掉末尾后能继续分解，但是在重构时会导致低频系数和高频系数长度不匹配的问题，此时就要确定分解级数j，在最开始将序列长度补为2^j的整数倍

注意2：而且增加序列长度的方式也很重要，直接补零会导致末尾有信号的突变，增添不必要的高频成分，所以采取镜像补充的方式。

4.音频信号测试

Daubechies(dbN)小波系

1.介绍

多贝西小波（英语：Daubechies Wavelet），是以比利时女性物理暨数学家英格丽·多贝西（Ingrid Daubechies）的名字命名之一种小波函数，当初英格丽·多贝西发现了一种具有阶层（hierarchy）性质的小波，便将此小波以她的名字命名。多贝西小波主要应用在离散型的小波转换，是最常使用到的小波转换，通常使用在数位信号分析、信号压缩跟噪声去除。

一般而言的离散小波转换通常是以正交小波（orthogonal wavelet）为基底，而多贝西小波也是一种正交小波。由于它很容易经由快速小波转换（fast wavelet transform（FWT））实现，所以常会放在数位信号处理的教科书中教学。

对于有限长度的小波，应用于快速小波转换（fast wavelet transform（FWT））时，会有两个实数组成的数列：一是作为高通滤波器的系数，称作小波滤波器（wavelet filter,也称为mother wavelet）；二是低通滤波器的系数，称作调整滤波器（scaling filter,也称为father wavelet）。

只要确定了滤波器系数，剩下的就直接与信号卷积即可。一般的滤波器是涉及频率响应再结合相位，用傅里叶逆变换回去。

小波变换的滤波器系数需要满足以下条件：

1.紧支撑性

2.小波的移动：两个两个向前移动，由下采样性质决定

3.正交性：利用正交性可以由母小波构造父小波，滤波器系数翻转然后隔一个数乘-1，来保证母小波和父小波是正交的

4.归一化：如果没有做归一化，得到的信号将越来越大

做了归一化之后，能够还原得到信号

5.零均值：简单理解为滤波器系数相加为0，这个系数只要求高通滤波器系数之和为0即可

6.消失矩：幂函数可以定义数据变化的剧烈程度，能抹平这种变化的能力强弱叫做函数的消失矩，如果与n次方函数点乘为0，那么消失矩为n+1。这里的n就是dbn系列里的n，如db4的含义是，这个小波的消失矩为4

7.离散小波的长度：正交对称性决定，离散小波长度应为偶数

8.如果要构造小波函数，可以按照以上条件，然后就变成了一个解方程组的问题，在系数慢慢增加后，就需要用到矩阵论中的逼近来求解。

得到了多贝西小波变换系数表

这里有一个疑问：为什么python运行得到的小波变换系数和表里的不一样，而且不满足零均值。

解答：零均值只需要高通滤波器满足即可，构造是以母小波为准。正交系数表里的数没有进行归一化，仅仅正交，归一化后就和python得到的一样了。归一化的操作是先求平方的和再开方，再用所有数除以这个得到的数就能归一化了。得到的与python运行出的系数是一样的。

2.分解代码

def db4_wavelet_decomposition(signal, signal_time=1, decomposition_level=1):"""对信号进行db4小波变化的分解:param f: 输入函数:param t: 信号时间:param level: 分解级数:return:approx, detail,dn"""# ==========1.设置滤波器系数==========# db4 低通滤波器 (8-tap)h = [-0.010597401785069032,0.0328830116668852,0.030841381835560764,-0.18703481171909309,-0.027983769416859854,0.6308807679298589,0.7148465705529157,0.2303778133088965]# db4 高通滤波器 (g = (-1)^k * h[7-k])g = [(-1) ** k * h[7 - k] for k in range(8)]n = len(signal)approx = np.zeros(n // 2)detail = np.zeros(n // 2)dn = [None] * decomposition_level  # 用于存储每一级高频系数#========== 开始逐级分解 ==========for level in range(decomposition_level):if level == 0:n = len(signal)print(n)fun = signaldn[0] = np.array([signal])else:n = len(approx)print(n)fun = approxdn[level] = np.array([detail])print(f"\n=== 正在进行第 {level + 1} 级分解 ===")print(f"当前信号长度: {n}")# 保存an之后创建新的an，dn用于存放下一级approx = np.zeros(n // 2)detail = np.zeros(n // 2)# 判断是否为偶数序列if n % 2 == 0:n = nelse:if n < 2:print("函数无法分解，程序终止")return False  # 返回 False 并跳出函数else:n = n - 1print('开始对序列进行db4小波变换，将其分解为低频信号与高频信号')# 先卷积再下采样for i in range(n):y = np.zeros(8)a = 0  # 近似系数累加器d = 0  # 细节系数累加器for m in range(8):if (i + m) >= n:y[m] = 0else:y[m] = fun[i + m]a += y[m] * h[m]d += y[m] * g[m]if i % 2 == 0:approx[i // 2] += adetail[i // 2] += d# 保留每一级的细节系数dn[level] = detailplt.subplot(3, 1, 1)plt.plot(np.linspace(0, signal_time, len(signal)), signal, label=f'Original Signal')plt.legend()plt.subplot(3, 1, 2)plt.plot(np.linspace(0, signal_time, len(approx)), approx, label=f'Approximation ')plt.legend()plt.subplot(3, 1, 3)plt.plot(np.linspace(0, signal_time, len(detail)), detail, label=f'detail ')plt.legend()plt.tight_layout()plt.show()return approx, detail,dn

如图是对信号进行了三级分解后的图像：

分析原信号的频率成分和三级分解后低频部分的频率成分：

代码如下：

def gffly(signal: np.ndarray,duration: float,plot: bool = True,title: Optional[str] = None
):"""简化的傅里叶变换频谱分析函数参数:----------signal : np.ndarray输入的一维时序信号duration : float信号时长(秒)plot : bool, optional是否绘制时域和频域图形(默认True)title : str, optional图形标题(默认None)返回:-------frequencies : np.ndarray频率数组(Hz)magnitude : np.ndarray幅度谱数组peaks : np.ndarray检测到的峰值频率数组(Hz)"""n_samples = len(signal)sample_rate = n_samples / duration  # 计算采样率t = np.linspace(0, duration, n_samples, endpoint=False)  # 时间轴# 应用汉宁窗(减少频谱泄漏)window = np.hanning(n_samples)signal_windowed = signal * window# 计算FFTfft_result = np.fft.fft(signal_windowed)frequencies = np.fft.fftfreq(n_samples, 1 / sample_rate)# 计算幅度谱magnitude = np.abs(fft_result)# 只取正频率部分half_n = n_samples // 2frequencies_half = frequencies[:half_n]magnitude_half = magnitude[:half_n]# 自动峰值检测(幅度大于最大幅度20%的峰值)peak_threshold = 0.2 * np.max(magnitude_half)peak_indices = np.where(magnitude_half > peak_threshold)[0]# 筛选真正的峰值(比相邻点高)true_peaks = []for i in peak_indices:if i == 0 or i == len(magnitude_half) - 1:continueif magnitude_half[i] > magnitude_half[i - 1] and magnitude_half[i] > magnitude_half[i + 1]:true_peaks.append(i)peaks = frequencies_half[true_peaks]# 绘图if plot:plt.figure(figsize=(12, 8))# 时域信号plt.subplot(2, 1, 1)plt.plot(t, signal)plt.title('Time Domain Signal' if title is None else f'Time Domain - {title}')plt.xlabel('Time [s]')plt.ylabel('Amplitude')plt.grid(True)# 频域幅度谱plt.subplot(2, 1, 2)plt.plot(frequencies_half, magnitude_half)plt.title('Frequency Spectrum' if title is None else f'Frequency Spectrum - {title}')plt.xlabel('Frequency [Hz]')plt.ylabel('Magnitude')plt.grid(True)# 标记检测到的峰值for peak in peaks:plt.annotate(f'{peak:.1f} Hz',xy=(peak, magnitude_half[np.where(frequencies_half == peak)[0][0]]),xytext=(10, 10), textcoords='offset points',arrowprops=dict(arrowstyle='->'))plt.tight_layout()plt.show()return frequencies_half, magnitude_half, peaks

跑出来的结果，左边是原信号的频率成分分析，右边是三级滤波之后的低频部分，可以看到低频部分主要是3Hz的成分，100Hz成分只有非常少的部分了，所以起到了滤波的作用

原信号的频率成分分析	三级滤波之后的低频部分

而高频部分产生的是25Hz的频率成分

3.重构代码

采用分解后的低频部分和高频部分对信号进行重构，重构滤波器就是将分解滤波器顺序反过来。

由这一级的两个系数重构得到上一级的低频系数

def db4_wavelet_reconstruction(approx, detail, duration=1):"""对db4小波变换进行重构:param approx::param detail::param duration::return:approx1"""h = [0.2303778133088965,0.7148465705529157,0.6308807679298589,-0.027983769416859854,-0.18703481171909309,0.030841381835560764,0.0328830116668852,-0.010597401785069032]# db4 高通滤波器 (g = (-1)^k * h[7-k])g = [(-1) ** k * h[7 - k] for k in range(8)]# 检查输入长度是否一致if len(approx) != len(detail):raise ValueError("近似系数和细节系数长度必须相同")n = len(approx)ann = np.zeros(2 * n)dnn = np.zeros(2 * n)approx1 = np.zeros(2 * n)# 对an和dn补零for i in range(n):ann[2 * i] = an[i]dnn[2 * i] = dn[i]print('开始进行小波变换重构')for i in range(2 * n):y = np.zeros(8)z = np.zeros(8)for m in range(8):if (i + m) >= 2 * n:y[m] = 0z[m] = 0else:y[m] = ann[i + m]z[m] = dnn[i + m]a = y * hd = z * gfor m in range(8):approx1[i] += a[m] + d[m]plt.subplot(3, 1, 1)plt.plot(np.linspace(0, duration, len(approx1)), approx1, label=f'Original Signal')plt.legend()plt.subplot(3, 1, 2)plt.plot(np.linspace(0, duration, len(approx)), approx, label=f'Approximation ')plt.legend()plt.subplot(3, 1, 3)plt.plot(np.linspace(0, duration, len(detail)), detail, label=f'detail ')plt.legend()plt.tight_layout()plt.show()return approx1

得到的重构图像如下：

这是没有进行滤波的重构，还可以增添阈值，对高频部分进行处理，这样重构的信号会有滤波的效果。

4.优化代码

3.真实信号测试

测量的加速度计和陀螺仪的真实信号，或者是音频信号

1.信号读取

需要处理的信号是以序列的形式存放，若信号为txt格式，第一列代表时间轴，第二列代表数据大小，信号时间轴以毫秒为单位读取数据，需要转换成序列信号，同时获取信号的时间作为返回值，便于下一步处理，wav为音频格式，两种信号是这样的：

txt	wav

读取txt格式信号：

def read_txt(file_path):"""读取不规则格式的TXT文件，第一列为毫秒时间戳参数:file_path: TXT文件路径返回:tuple: (signal_array, time_seconds_array, total_duration_seconds)- signal_array: 信号值数组- time_seconds_array: 时间数组(秒为单位)- total_duration_seconds: 信号总时长(秒)"""time_ms = []  # 存储毫秒时间戳signal = []  # 存储信号值with open(file_path, 'r') as file:for line in file:line = line.strip()# 跳过空行和注释行（以#开头的行）if not line or line.startswith('#'):continue# 分割行内容parts = line.split()if len(parts) >= 2:try:# 第一列为时间戳(毫秒)，最后一列为信号值timestamp = float(parts[0])value = float(parts[-1])time_ms.append(timestamp)signal.append(value)except ValueError:continue# 转换为numpy数组time_ms = np.array(time_ms)signal = np.array(signal)# 将毫秒转换为秒time_seconds = time_ms / 1000.0# 计算总时长(秒)if len(time_seconds) > 1:total_duration = time_seconds[-1] - time_seconds[0]else:total_duration = 0.0return signal, time_seconds, total_duration

读取wav格式信号：

    wave_path = r'D:/project/signalProcessing/testdata/cay.wav'file = wave.open(wave_path)# for item in enumerate(WAVE.getparams()):#     print(item)a = file.getparams().nframes  # 帧总数f = file.getparams().framerate  # 采样频率sample_frequency, audio_sequence = wavfile.read(wave_path)sample_time = 1 / f  # 采样点的时间间隔time = a / f  # 声音信号的长度# 读取原始音频（保留原始数据格式）data, original_sr = sf.read(wave_path)# print(audio_sequence)  # 声音信号每一帧的“大小”x_seq = np.arange(0, time, sample_time)

这里audio_sequence就是读取到的信号序列，sample_time就是信号时间长度。

2.信号处理（分解，滤波，重构）

3.信号分析对比（频谱分析，信噪比等）

1.信号相似性度量指标

最简单的，可以直接作差，观察重构信号与原信号，同时找出二者差值的最大值，代码如下

def compare(signal1,signal2,t = 1):"""输入两个信号，一个是处理之前的信号，另一个是滤波之后的信号将二者进行相似度对比:param signal1::param signal2::return:dif"""n = len(signal1)if len(signal1) != len(signal2):print('信号一长度',n,'信号二长度',len(signal2))raise ValueError("输入两信号长度必须相同")dif = np.zeros(n)# ========== 作差 ===========for i in range(n):dif[i] = signal1[i] - signal2[i]print('最大差值为:',max(abs(dif)))# ========== 可视化 ==========plt.subplot(3, 1, 1)plt.plot(np.linspace(0, t, n), signal1, label=f'signal1')plt.legend()plt.subplot(3, 1, 2)plt.plot(np.linspace(0, t, n), signal2, label=f'signal2')plt.legend()plt.subplot(3, 1, 3)plt.plot(np.linspace(0, t, n), dif, label=f'dif')plt.legend()plt.tight_layout()plt.show()return dif

用正弦波信号测试