信号各种变换 python实现 DFT,STFT,CWT,DWT

数据读取

傅里叶变换

短时傅里叶变换

连续小波变换

离散小波变换

import numpy as np
import matplotlib.pyplot as plt
import pywt.data
from scipy.fftpack import fft,ifft
from scipy import signal

数据读取

首先我们有这么一个信号，下图是数据格式，即一列一个数：

我们读取以后只要前500个数：

a = []
b = []
with open("ECG1_1.txt") as file_obj:for content in file_obj:a.append(int(content))
for i in range(500):b.append(a[i])
plt.plot(b)
plt.show()

傅里叶变换

fft_b=fft(b)abs_b = np.abs(fft_b)  # 取复数的模
angle_b = np.angle(fft_b)  # 取复数的幅角plt.figure()
plt.plot(abs_b)
plt.figure()
plt.plot(angle_b)plt.show()

其中，上图是频率振幅，下图是频率幅角：

短时傅里叶变换

#  fs:时间序列的采样频率,  nperseg:每个段的长度   noverlap:段之间重叠的点数。如果没有则noverlap=nperseg/2
f, t, nd = signal.stft(b ,fs = 1.0,window ='hann',nperseg = 150,noverlap = 50)
plt.pcolormesh(t, f, np.abs(nd), vmin = 0, vmax = 4)
plt.title('STFT')
plt.ylabel('frequency')
plt.xlabel('time')
plt.show()

显示结果如下：

但是这个着实不太好分析，我们换个简单一点的函数：

aa = []
for i in range(200):aa.append(np.sin(0.3*np.pi*i))
for i in range(200):aa.append(np.sin(0.13*np.pi*i))
for i in range(200):aa.append(np.sin(0.05*np.pi*i))
plt.plot(aa)
plt.show()#  fs:时间序列的采样频率,  nperseg:每个段的长度   noverlap:段之间重叠的点数。如果没有则noverlap=nperseg/2
f, t, nd = signal.stft(aa ,fs = 1.0,window ='hann',nperseg = 150,noverlap = 50)
plt.pcolormesh(t, f, np.abs(nd), vmin = 0, vmax = 4)
plt.title('STFT')
plt.ylabel('frequency')
plt.xlabel('time')
plt.show()

显示为：

频率显示也非常正常：一开始频率很高，然后第二段降低，第三段最低，同时重叠区域有两段频率：

设置窗不重叠，则显示为：

连续小波变换

sampling_rate = 1024
wavename = 'cgau8'
totalscal = 256
# 中心频率
fc = pywt.central_frequency(wavename)
# 计算对应频率的小波尺度
cparam = 2 * fc * totalscal
scales = cparam / np.arange(totalscal, 0, -1)
[cwtmatr, frequencies] = pywt.cwt(aa, scales, wavename, 1.0 / sampling_rate)
plt.figure(figsize=(8, 4))
plt.subplot(211)
t = np.arange(0, 600, 1.0)
plt.plot(t, aa)
plt.xlabel(u"time(s)")
plt.subplot(212)
plt.contourf(t, frequencies, abs(cwtmatr))
plt.ylabel(u"freq(Hz)")
plt.xlabel(u"time(s)")
plt.subplots_adjust(hspace=0.4)
plt.show()

具体的解释可以参考我的网站中：

小波分析 – Dezeming Familyhttps://dezeming.top/?page_id=1019的关于Python小波分析的专栏。

修改一下原始信号，变成中间高频，两边低频，得到结果：

离散小波变换

wavename = 'db5'
cA, cD = pywt.dwt(aa, wavename)
ya = pywt.idwt(cA, None, wavename,'smooth') # approximated component
yd = pywt.idwt(None, cD, wavename,'smooth') # detailed component
x = range(len(aa))
plt.figure(figsize=(12,9))
plt.subplot(311)
plt.plot(x, aa)
plt.title('original signal')
plt.subplot(312)
plt.plot(x, ya)
plt.title('approximated component')
plt.subplot(313)
plt.plot(x, yd)
plt.title('detailed component')
plt.tight_layout()
plt.show()

离散小波变换DWT：

CWT的“连续性”，以及它与离散小波变换的区别，是它运行的尺度和位置集。与离散小波变换不同，CWT可以在每一个尺度上进行操作，从原始信号的尺度到某个最大尺度。（当然对于计算机来说，也是从中抽取一定数量的离散尺度）

而且CWT在移位方面也是连续的，即在计算过程中，分析函数的整个域上平滑地移位（对于计算机来说，也是根据时间分辨率来离散的移位）。

上面的函数把信号分成了低频近似和高频细节两个部分：