👨🏼‍🎓 🛀🏻 🥠 小波分析。第二部分 🥟 👩🏽‍🤝‍👨🏻 🎯

引言

该出版物讨论了时间序列的小波分析。小波变换的主要思想与许多时间序列的细节相对应，这些细节证明了其主要特征随时间的演变-谐波分量的平均值，方差，周期，幅度和相位。在各个知识领域研究的绝大多数过程都具有上述特征。

该出版物的目的是描述使用PyWavelets库对时间序列进行连续小波变换的方法。

一点历史

二十世纪70年代后期的地球物理学家D. Morlaix。我遇到了分析来自地震传感器的信号的问题，这些信号在短时间内包含高频成分（地震活动），在很长一段时间内包含低频成分（地壳的平静状态）。窗口傅立叶变换允许您分析高频分量或低频分量，但不能一次分析两个分量。

因此，提出了一种分析方法，其中低频的窗函数的宽度增加，而高频的窗函数的宽度减小。作为一个生成函数（所谓的缩放函数-缩放函数，scalet）的拉伸（压缩）和时间偏移的结果，获得了新的窗口变换。 D. Morlaix将此生成函数称为小波。

小波D.Morlaix

from pylab import* import scaleogram as scg axes = scg.plot_wav('cmor1-1.5', figsize=(14,3)) show()

连续小波变换（CWT）

单级小波变换：

系数，频率= pywt.cwt（数据，标度，小波）

其中：

数据： （array_like） -输入信号;

scales （array_like） ：-要使用的小波尺度。您可以使用比率f = scale2frequency（比例，小波）/ sample_period来确定物理频率f是多少。此处，f以秒为单位时的samples_period为赫兹；

小波 （小波对象或名称） ：-要使用的小波；

采样周期（浮点）：-输出频率的采样周期（可选参数）。为系数计算的值与sample_period无关（即，比例未通过周期采样进行缩放）。

coefs （array_like） ：-连续小波-给定尺度和小波的输入信号变换；

频率 （array_like） ：-如果采样周期的单位是秒并已设置，则频率以赫兹显示。否则，假定采样周期为1。

注意：系数数组的大小取决于输入数组的长度和给定比例的长度。

获取可用的cwt兼容小波名称列表：

 >>> import pywt >>> wavlist = pywt.wavelist(kind='continuous') >>> wavlist

['cgau1'，'cgau2'，'cgau3'，'cgau4'，'cgau5'，'cgau6'，'cgau7'，'cgau8'，'cmor'，'fbsp'，'gaus1'，'gaus2'，' gaus3'，'gaus4'，'gaus5'，'gaus6'，'gaus7'，'gaus8'，'mexh'，'morl'，'shan']]

对于小波分析，选择母小波是关键任务之一。因此，在每次选择小波之前，只需使用以下程序即可确定小波的动态属性：

清单1

 import numpy as np import pywt import matplotlib.pyplot as plt vav='cmor1.5-1.0' wav = pywt.ContinuousWavelet(vav) #  ,      print("       [{}, {}]".format( wav.lower_bound, wav.upper_bound)) width = wav.upper_bound - wav.lower_bound scales = [1, 2, 3, 4, 10, 15] max_len = int(np.max(scales)*width + 1) t = np.arange(max_len) fig, axes = plt.subplots(len(scales), 2, figsize=(12, 6)) for n, scale in enumerate(scales): #       cwt int_psi, x = pywt.integrate_wavelet(wav, precision=10) step = x[1] - x[0] j = np.floor( np.arange(scale * width + 1) / (scale * step)) if np.max(j) >= np.size(int_psi): j = np.delete(j, np.where((j >= np.size(int_psi)))[0]) j = j.astype(np.int) # normalize int_psi     int_psi /= np.abs(int_psi).max() #     filt = int_psi[j][::-1] # CWT          #          . nt = len(filt) t = np.linspace(-nt//2, nt//2, nt) axes[n, 0].plot(t, filt.real, t, filt.imag) axes[n, 0].set_xlim([-max_len//2, max_len//2]) axes[n, 0].set_ylim([-1, 1]) axes[n, 0].text(50, 0.35, 'scale = {}'.format(scale)) f = np.linspace(-np.pi, np.pi, max_len) filt_fft = np.fft.fftshift(np.fft.fft(filt, n=max_len)) filt_fft /= np.abs(filt_fft).max() axes[n, 1].plot(f, np.abs(filt_fft)**2) axes[n, 1].set_xlim([-np.pi, np.pi]) axes[n, 1].set_ylim([0, 1]) axes[n, 1].set_xticks([-np.pi, 0, np.pi]) axes[n, 1].set_xticklabels([r'$-\pi$', '0', r'$\pi$']) axes[n, 1].grid(True, axis='x') axes[n, 1].text(np.pi/2, 0.5, 'scale = {}'.format(scale)) axes[n, 0].set_xlabel('time (samples)') axes[n, 1].set_xlabel('frequency (radians)') axes[0, 0].legend(['real', 'imaginary'], loc='upper left') axes[0, 1].legend(['Power'], loc='upper left') axes[0, 0].set_title('filter: wavelet - %s'%vav) axes[0, 1].set_title(r'|FFT(filter)|$^2$') plt.show()

接下来，我们将使用上述程序考虑小波函数及其属性：

墨西哥帽小波混合：

$\ varphi \左（t \右）= \ frac {2} {\ sqrt {3} \ sqrt [4] {\ pi}} \ cdot exp ^ {-\ frac {t ^ {2}} {2} } \ cdot \ left（1-t ^ {2} \ right）$

小波“ morl” Morlaix：

$\ varphi \ left（t \ right）= exp ^ {-\ frac {t ^ {2}} {2}} \ cdot cos（5t）$

复Morlet小波“ cmorB-C”

$\ varphi \ left（t \ right）= \ frac {1} {\ sqrt {\ pi} B} exp ^ {-\ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

其中B是带宽； C是中心频率。

在下文中（无特殊说明），B，C的值设置为浮点数。

高斯小波“ gausP”

$\ varphi \左（t \右）= C \ cdot exp ^ {-t ^ {2}}$

其中：P-小波名称中插入1到8的整数，例如“ gaus8”； C-内置标准化常数。

集成高斯小波“ cgauP”

$\ varphi \ left（t \ right）= C \ cdot exp ^ {-t ^ {2}} \ cdot exp ^ {-jt}$

其中：P-小波名称中插入1到8的整数，例如“gaus8”，它对应于小波函数的导数的阶数； C-内置标准化常数。

香农小波“ shanB-C”

$\ varphi \ left（t \ right）= \ sqrt {B} \ cdot \ frac {sin（\ pi B t）} {\ pi B T} \ cdot exp ^ {j 2 pi C t}$

其中：B是带宽； C是中心频率。

PyWavelets中的CWT适用于离散数据-具有小波积分样本的卷积。如果比例尺太小，那么我们将得到一个采样不足的离散滤波器，从而导致平滑，如小波“ cmor1.5-1.0”的图形所示。

在左列中，图表显示了在卷积中使用的各种比例的离散滤波器。右列是每个滤波器的相应傅立叶功率谱。对于图“ cmor1.5-1.0”，使用比例1和2可以看出，由于违反了奈奎斯特约束，因此发生了平滑。当选择C和B时，指示的现象可能在其他小波中发生，因此，在使用小波时，必须使用程序-清单1。

要将给定的标度与特定的信号频率相关联，必须知道信号采样周期。使用pywt.scale2frequency（）函数，可以将刻度列表转换为相应的频率。比例尺的正确选择取决于所选的小波，因此应使用pywt.scale2frequency（）来了解信号的相应频率范围。

 >>> import pywt >>> dt = 0.01 # 100 Hz sampling >>> frequencies = pywt.scale2frequency('cmor1.5-1.0', [1, 2, 3, 4]) / dt >>> frequencies

数组（[100.，50.，33.33333333，25.]）

小波-使用CWT PyWavelet进行时间序列分析

el-Nino数据集是用于跟踪El Nino的时间序列数据集，包含从1871年到1997年的每季度海表温度测量值。为了了解比例尺图表的功能，让我们将其可视化为el-Nino数据集以及原始时间序列数据及其傅里叶变换。

对于时间序列的小波分析，必须执行以下步骤：

1.选择母小波：选择复数Morlet小波“ cmorB-C”：

$\ varphi \ left（t \ right）= \ frac {1} {\ sqrt {\ pi} B} exp ^ {-\ frac {t ^ {2}} {B}} \ cdot exp ^ {j2 \ pi C t}$

带宽-B和中心频率-下一步将确定带宽。

2.通过使我们的时间序列的dt = 0.25来确定中心频率：

 import pywt dt = 0.25 # 4 Hz sampling scale=range(1,4) frequencies = pywt.scale2frequency('cmor1.0-0.5', scale) / dt print(frequencies)

[2。 1. 0.66666667]

3.我们找到了母波cmor1.0-0.5的傅立叶变换（我们使用清单1）：
连续小波将在整个范围内进行评估[-8.0，8.0]

4.选择时间序列的数据：

 pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|")

数据是1871年至1997年的季度海表温度测量值。对于此数据：t0 = 1871 dt = 0.25

5.我们在一个图表上建立一个信号及其移动平均值的时间序列：

清单2

 from numpy import * from scipy import * import pandas as pd from pylab import * import pywt def get_ave_values(xvalues, yvalues, n = 5): signal_length = len(xvalues) if signal_length % n == 0: padding_length = 0 else: padding_length = n - signal_length//n % n xarr = array(xvalues) yarr = array(yvalues) xarr.resize(signal_length//n, n) yarr.resize(signal_length//n, n) xarr_reshaped = xarr.reshape((-1,n)) yarr_reshaped = yarr.reshape((-1,n)) x_ave = xarr_reshaped[:,0] y_ave = nanmean(yarr_reshaped, axis=1) return x_ave, y_ave def plot_signal_plus_average(time, signal, average_over = 5): fig, ax = subplots(figsize=(15, 3)) time_ave, signal_ave = get_ave_values(time, signal, average_over) ax.plot(time, signal, label='') ax.plot(time_ave, signal_ave, label = '   (n={})'.format(5)) ax.set_xlim([time[0], time[-1]]) ax.set_ylabel(' ', fontsize=18) ax.set_title(' +   ', fontsize=18) ax.set_xlabel('', fontsize=18) ax.legend() show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_signal_plus_average(time, signal)

6.我们从时间序列中进行傅立叶变换和频谱模态：

清单3

 from numpy import * from scipy import * import pandas as pd from pylab import * import pywt def get_ave_values(xvalues, yvalues, n = 5): signal_length = len(xvalues) if signal_length % n == 0: padding_length = 0 else: padding_length = n - signal_length//n % n xarr = array(xvalues) yarr = array(yvalues) xarr.resize(signal_length//n, n) yarr.resize(signal_length//n, n) xarr_reshaped = xarr.reshape((-1,n)) yarr_reshaped = yarr.reshape((-1,n)) x_ave = xarr_reshaped[:,0] y_ave = nanmean(yarr_reshaped, axis=1) return x_ave, y_ave def get_fft_values(y_values, T, N, f_s): f_values = linspace(0.0, 1.0/(2.0*T), N//2) fft_values_ = fft(y_values) fft_values = 2.0/N * abs(fft_values_[0:N//2]) return f_values, fft_values def plot_fft_plus_power(time, signal): dt = time[1] - time[0] N = len(signal) fs = 1/dt fig, ax = subplots(figsize=(15, 3)) variance = std(signal)**2 f_values, fft_values = get_fft_values(signal, dt, N, fs) fft_power = variance * abs(fft_values) ** 2 # FFT power spectrum ax.plot(f_values, fft_values, 'r-', label='FFT ') ax.plot(f_values, fft_power, 'k--', linewidth=1, label=' ') ax.set_xlabel('[ / ]', fontsize=18) ax.set_ylabel('', fontsize=18) ax.legend() show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_fft_plus_power(time, signal)

7.我们确定尺度：尺度=范围（1，128）; 级别= [2 **-4，4 **-3，2 **-2，2 **-1，2 ** 0，2 ** 1，2 ** 2，2 ** 3]。

8.建立一个比例表：

清单4

 from numpy import * import pandas as pd from pylab import * import pywt def plot_wavelet(time, signal, scales, waveletname = 'cmor1.0-0.4', cmap = plt.cm.seismic, title = '-( ) ', ylabel = ' ()', xlabel = ''): dt = time[1] - time[0] [coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt) power = (abs(coefficients)) ** 2 period = 1. / frequencies levels = [2**-4 , 2**-3 , 2**-2 , 2**-1 , 2**0 , 2**1 , 2**2 , 2**3] contourlevels = log2(levels) fig, ax = subplots(figsize=(15, 10)) im = ax.contourf(time, log2(period), log2(power), contourlevels, extend='both',cmap=cmap) ax.set_title(title, fontsize=20) ax.set_ylabel(ylabel, fontsize=18) ax.set_xlabel(xlabel, fontsize=18) yticks = 2**arange(np.ceil(log2(period.min())), ceil(log2(period.max()))) ax.set_yticks(log2(yticks)) ax.set_yticklabels(yticks) ax.invert_yaxis() ylim = ax.get_ylim() ax.set_ylim(ylim[0], -1) cbar_ax = fig.add_axes([0.95, 0.5, 0.03, 0.25]) fig.colorbar(im, cax=cbar_ax, orientation="vertical") show() df_nino = pd.read_csv("http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat", sep = "|") N = df_nino.shape[0] t0=1871 dt=0.25 time = arange(0, N) * dt + t0 signal = df_nino.values.squeeze() scales = arange(1, 128) plot_wavelet(time, signal, scales)

比例表显示，大多数频谱功率集中在2-8年内，这对应于0.125-0.5 Hz的频率。在傅立叶频谱中，时尚也集中在这些频率值上。但是，小波变换也提供了临时信息，而傅立叶变换则没有。

例如，在比例图上，我们看到1920年之前波动很大，而1960年至1990年之间波动很小。我们还可以看到，随着时间的流逝，从较短的时期向较长的时期转变。

使用比例尺模块

Scaleogram是基于PyWavelets库通过连续小波变换分析一维数据的便捷工具。该模块旨在为快速数据分析提供可靠的工具。

比例图具有以下功能：

简单的初学者通话签名
可读轴和纯matplotlib集成
缩放，频谱过滤器，颜色条实现等许多选项...
根据标记支持频率和频率单位
速度，使用N * log（N）算法进行转换
可移植性：已通过python2.7和python3.7测试
详细的错误消息和文档以及示例

但是，在使用示例中，数据访问未正确记录：

 nino3 = "http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat"

出现以下警告：

nino3 = pd.read_table（nino3）
导致警告：
警告（来自警告模块）：
文件“ C：/Users/User/Desktop/wavelet/pr.py”，第7行
nino3 = pd.read_table（nino3）
FutureWarning：不建议使用read_table，请改用read_csv，并传递sep ='\ t'

数据访问应这样写：

 url="http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|")

为了显示比例尺的使用并将结果与上一个示例的结果进行比较，我们使用相同的数据-对1871年至1997年的海面温度进行季度测量。对于此数据：t0 = 1871 dt = 0.25

清单4

 import pandas as pd import pywt from numpy import* import scaleogram as scg from pylab import* url="sst_nino3.dat" nino3 = pd.read_csv(url, sep = "|") data = nino3.values.squeeze() N = data.size t0 = 1871; dt = 0.25 time = t0 + arange(len(data))*dt wavelet = 'cmor1-0.5' ax = scg.cws(time, data, scales=arange(1, 128), wavelet=wavelet, figsize=(14, 3), cmap="jet", cbar=None, ylabel=' ()', xlabel=" []", yscale="log", title='-  \n( )') ticks = ax.set_yticks([2,4,8, 16,32]) ticks = ax.set_yticklabels([2,4,8, 16,32]) show()

如果我们将比例尺图表与获得的扫描结果进行比较，则除了调色板外，它们是相同的，因此在时间序列中显示出相同的动态。
清单4比清单3更简单，并且具有更好的性能。

当根据傅立叶的频谱和功率谱不允许获得有关时间序列动态的附加信息时，则最好使用比例图。

结论

给出了通过PyWavelets库对时间序列进行小波分析的示例。

小波分析。 第二部分

引言

连续小波变换（CWT）

小波-使用CWT PyWavelet进行时间序列分析

使用比例尺模块

结论

More articles:

小波分析。第二部分