从‘厄尔尼诺’数据到你的数据用PyWavelets绘制专业级CWT时频图避坑指南在时间序列分析领域小波变换正逐渐成为研究非平稳信号的利器。与传统的傅里叶变换相比连续小波变换(CWT)能够同时捕捉信号的时域和频域特征特别适合分析那些频率随时间变化的复杂信号。本文将带你深入PyWavelets库的CWT实现从厄尔尼诺气象数据入手最终应用到你的股票价格、设备振动等实际数据中。1. 理解小波变换的核心概念小波变换之所以强大源于其独特的时频局部化能力。想象一下分析一段音乐信号傅里叶变换能告诉你有哪些音符但不知道它们何时出现而小波变换不仅能识别音符还能精确定位每个音符的演奏时刻。小波与正弦波的关键区别正弦波从负无穷延伸到正无穷没有时间局部性小波是有限长度的波形具有明确的开始和结束小波变换通过平移和缩放操作实现了多分辨率分析PyWavelets库中pywt.cwt()函数的数学本质是卷积运算CWT(a,τ) ∫ f(t)ψ*((t-τ)/a) dt其中a是尺度参数τ是平移参数ψ是小波函数。提示尺度与频率成反比关系较大尺度对应较低频率较小尺度对应较高频率2. 配置CWT分析的关键参数2.1 母小波的选择艺术PyWavelets提供了14种小波家族选择合适的小波对分析结果至关重要小波类型特点适用场景Morlet (morl)对称、非正交时频分析通用选择Mexican Hat (mexh)二阶导数形式边缘检测Daubechies (db)紧支撑正交信号压缩Complex Morlet (cmor)复数小波相位分析# 查看可用连续小波 print(pywt.wavelist(kindcontinuous))2.2 尺度参数的智能生成尺度参数决定了分析的频率范围常见错误是直接使用线性间隔的尺度import numpy as np def generate_scales(wavelet, sr, num_scales150): fc pywt.central_frequency(wavelet) # 中心频率 cparam 2 * fc * num_scales scales cparam / np.arange(num_scales, 1, -1) return scales # 示例为采样率128Hz的信号生成尺度 scales generate_scales(morl, sr128)2.3 采样率与频率转换频率轴标注不准确是常见问题正确的频率转换方法frequencies pywt.scale2frequency(morl, scales) * sr3. 绘制专业级时频图的实战技巧3.1 避免图像模糊的配置plt.figure(figsize(12, 6)) plt.contourf(t, frequencies, np.abs(cwtmatr), levels100, cmapjet, vmaxabs(cwtmatr).max(), vmin-abs(cwtmatr).max()) plt.colorbar() plt.yscale(log) # 对数频率轴更直观常见问题解决方案图像模糊 → 增加contourf的levels参数频率轴不直观 → 使用对数坐标颜色对比不足 → 调整vmin/vmax参数3.2 3D时频图的高级呈现from mpl_toolkits.mplot3d import Axes3D fig plt.figure(figsize(14, 8)) ax fig.add_subplot(111, projection3d) X, Y np.meshgrid(t, frequencies) ax.plot_surface(X, Y, np.abs(cwtmatr), cmapviridis, rstride1, cstride1) ax.view_init(30, 45)3.3 厄尔尼诺数据分析案例# 加载厄尔尼诺数据集 url http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat data pd.read_table(url, headerNone).values.flatten() # 设置时间轴 years np.arange(1871, 1997.25, 0.25) # 执行CWT分析 scales np.arange(1, 128) coefficients, freqs pywt.cwt(data, scales, morl, sampling_period0.25) # 可视化 plt.figure(figsize(15, 5)) plt.contourf(years, freqs, np.abs(coefficients), levels100, cmapRdYlBu_r) plt.colorbar(labelMagnitude) plt.title(ENSO Index Wavelet Power Spectrum) plt.xlabel(Year) plt.ylabel(Period (years)) plt.yscale(log)4. 应用到你的数据完整工作流4.1 股票价格分析实战# 获取股票数据 import yfinance as yf msft yf.Ticker(MSFT) hist msft.history(period5y) # 预处理 close_prices hist[Close].values dates hist.index.to_pydatetime() # 执行CWT scales generate_scales(cmor1.5-1.0, sr252) # 252个交易日 coefficients, freqs pywt.cwt(close_prices, scales, cmor1.5-1.0) # 可视化 plt.figure(figsize(15, 5)) plt.contourf(dates, freqs, np.abs(coefficients), levels50, cmapviridis) plt.colorbar() plt.title(Microsoft Stock Price Wavelet Analysis) plt.ylabel(Frequency (1/year)) plt.yscale(log)4.2 工业设备振动监测# 模拟振动传感器数据 t np.linspace(0, 10, 5000) vibration np.sin(2*np.pi*50*t) * (t2) np.sin(2*np.pi*120*t) * (t5) # 故障特征分析 scales generate_scales(mexh, sr500) cwtmatr, freqs pywt.cwt(vibration, scales, mexh) # 专业可视化 fig, (ax1, ax2) plt.subplots(2, 1, figsize(12, 8)) ax1.plot(t, vibration) ax1.set_title(Raw Vibration Signal) img ax2.imshow(np.abs(cwtmatr), extent[t.min(), t.max(), freqs.min(), freqs.max()], cmaphot, aspectauto, vmaxabs(cwtmatr).max(), vmin0) fig.colorbar(img, axax2) ax2.set_title(Continuous Wavelet Transform) ax2.set_ylabel(Frequency [Hz])在工业应用中小波变换能清晰捕捉到设备故障发生时频率成分的变化比传统的FFT分析更能提前预警异常情况。
从‘厄尔尼诺’数据到你的数据：用PyWavelets绘制专业级CWT时频图避坑指南
从‘厄尔尼诺’数据到你的数据用PyWavelets绘制专业级CWT时频图避坑指南在时间序列分析领域小波变换正逐渐成为研究非平稳信号的利器。与传统的傅里叶变换相比连续小波变换(CWT)能够同时捕捉信号的时域和频域特征特别适合分析那些频率随时间变化的复杂信号。本文将带你深入PyWavelets库的CWT实现从厄尔尼诺气象数据入手最终应用到你的股票价格、设备振动等实际数据中。1. 理解小波变换的核心概念小波变换之所以强大源于其独特的时频局部化能力。想象一下分析一段音乐信号傅里叶变换能告诉你有哪些音符但不知道它们何时出现而小波变换不仅能识别音符还能精确定位每个音符的演奏时刻。小波与正弦波的关键区别正弦波从负无穷延伸到正无穷没有时间局部性小波是有限长度的波形具有明确的开始和结束小波变换通过平移和缩放操作实现了多分辨率分析PyWavelets库中pywt.cwt()函数的数学本质是卷积运算CWT(a,τ) ∫ f(t)ψ*((t-τ)/a) dt其中a是尺度参数τ是平移参数ψ是小波函数。提示尺度与频率成反比关系较大尺度对应较低频率较小尺度对应较高频率2. 配置CWT分析的关键参数2.1 母小波的选择艺术PyWavelets提供了14种小波家族选择合适的小波对分析结果至关重要小波类型特点适用场景Morlet (morl)对称、非正交时频分析通用选择Mexican Hat (mexh)二阶导数形式边缘检测Daubechies (db)紧支撑正交信号压缩Complex Morlet (cmor)复数小波相位分析# 查看可用连续小波 print(pywt.wavelist(kindcontinuous))2.2 尺度参数的智能生成尺度参数决定了分析的频率范围常见错误是直接使用线性间隔的尺度import numpy as np def generate_scales(wavelet, sr, num_scales150): fc pywt.central_frequency(wavelet) # 中心频率 cparam 2 * fc * num_scales scales cparam / np.arange(num_scales, 1, -1) return scales # 示例为采样率128Hz的信号生成尺度 scales generate_scales(morl, sr128)2.3 采样率与频率转换频率轴标注不准确是常见问题正确的频率转换方法frequencies pywt.scale2frequency(morl, scales) * sr3. 绘制专业级时频图的实战技巧3.1 避免图像模糊的配置plt.figure(figsize(12, 6)) plt.contourf(t, frequencies, np.abs(cwtmatr), levels100, cmapjet, vmaxabs(cwtmatr).max(), vmin-abs(cwtmatr).max()) plt.colorbar() plt.yscale(log) # 对数频率轴更直观常见问题解决方案图像模糊 → 增加contourf的levels参数频率轴不直观 → 使用对数坐标颜色对比不足 → 调整vmin/vmax参数3.2 3D时频图的高级呈现from mpl_toolkits.mplot3d import Axes3D fig plt.figure(figsize(14, 8)) ax fig.add_subplot(111, projection3d) X, Y np.meshgrid(t, frequencies) ax.plot_surface(X, Y, np.abs(cwtmatr), cmapviridis, rstride1, cstride1) ax.view_init(30, 45)3.3 厄尔尼诺数据分析案例# 加载厄尔尼诺数据集 url http://paos.colorado.edu/research/wavelets/wave_idl/sst_nino3.dat data pd.read_table(url, headerNone).values.flatten() # 设置时间轴 years np.arange(1871, 1997.25, 0.25) # 执行CWT分析 scales np.arange(1, 128) coefficients, freqs pywt.cwt(data, scales, morl, sampling_period0.25) # 可视化 plt.figure(figsize(15, 5)) plt.contourf(years, freqs, np.abs(coefficients), levels100, cmapRdYlBu_r) plt.colorbar(labelMagnitude) plt.title(ENSO Index Wavelet Power Spectrum) plt.xlabel(Year) plt.ylabel(Period (years)) plt.yscale(log)4. 应用到你的数据完整工作流4.1 股票价格分析实战# 获取股票数据 import yfinance as yf msft yf.Ticker(MSFT) hist msft.history(period5y) # 预处理 close_prices hist[Close].values dates hist.index.to_pydatetime() # 执行CWT scales generate_scales(cmor1.5-1.0, sr252) # 252个交易日 coefficients, freqs pywt.cwt(close_prices, scales, cmor1.5-1.0) # 可视化 plt.figure(figsize(15, 5)) plt.contourf(dates, freqs, np.abs(coefficients), levels50, cmapviridis) plt.colorbar() plt.title(Microsoft Stock Price Wavelet Analysis) plt.ylabel(Frequency (1/year)) plt.yscale(log)4.2 工业设备振动监测# 模拟振动传感器数据 t np.linspace(0, 10, 5000) vibration np.sin(2*np.pi*50*t) * (t2) np.sin(2*np.pi*120*t) * (t5) # 故障特征分析 scales generate_scales(mexh, sr500) cwtmatr, freqs pywt.cwt(vibration, scales, mexh) # 专业可视化 fig, (ax1, ax2) plt.subplots(2, 1, figsize(12, 8)) ax1.plot(t, vibration) ax1.set_title(Raw Vibration Signal) img ax2.imshow(np.abs(cwtmatr), extent[t.min(), t.max(), freqs.min(), freqs.max()], cmaphot, aspectauto, vmaxabs(cwtmatr).max(), vmin0) fig.colorbar(img, axax2) ax2.set_title(Continuous Wavelet Transform) ax2.set_ylabel(Frequency [Hz])在工业应用中小波变换能清晰捕捉到设备故障发生时频率成分的变化比传统的FFT分析更能提前预警异常情况。
