代码共享 | 使用Python执行快速傅里叶变换和连续小波变换

import pywtimport numpy as npimport matplotlib.pyplot as pltimport pandas as pdplt.style.use('seaborn')
dataset = "./monsoon.txt"df_nino = pd.read_table(dataset,header=None)df_nino

def plot_signal(time, signal, average_over=5, figname=None):    fig, ax = plt.subplots(figsize=(15, 3))    # txt是一个纯文本，行全部取，列的话只取第二列，1：    # 这里的time相当于x轴，signal是需要绘制的曲线    ax.plot(time, signal[:,1:], label='signal',lw=4.0)    ax.set_xlim([time[0], time[-1]])    # 设置x轴的范围，第一个时间到最后一个时间点    ax.set_ylabel('Signal Amplitude', fontsize=16)    ax.set_title('Signal + Time Average', fontsize=16)    ax.set_xlabel('Time', fontsize=16)    ax.legend()    plt.show()#     if not figname:#         plt.savefig('signal_plot.png', dpi=300, bbox_inches='tight')#     else:#         plt.savefig(figname, dpi=300, bbox_inches='tight')#     plt.close('all')
plot_signal(time, signal) #plot and label the axis

def get_fft_values(y_values, T, N, f_s):    # 获取时间序列的FFT变换值，快速傅里叶变换    #     f_values = np.linspace(0.0, 1.0/(2.0*T), N//2)    fft_values_ = np.fft.fft(y_values)    fft_values = 2.0/N * np.abs(fft_values_[0:N//2])    # FFT序列的计算    return f_values, fft_values

def plot_fft_plus_power(time, signal, figname=None):    dt = time[1] - time[0]    N = len(signal)    fs = 1/dt    #     fig, ax = plt.subplots(2, 1, figsize=(15, 6), sharex=True)    # 获取方差    variance = np.std(signal[:,1


    
:])**2    f_values, fft_values = get_fft_values(signal[:,1:], dt, N, fs)    # 获取f_values之后计算FFT谱值序列=序列方差*fft值的绝对值在二次幂    fft_power = variance * abs(fft_values) ** 2  # FFT power spectrum    ax[0].plot(f_values, fft_values, 'r-', label='Fourier Transform')    ax[1].plot(f_values, fft_power, 'b-',               linewidth=1, label='FFT Power Spectrum')    ax[1].set_xlabel('Frequency [Hz / year]', fontsize=18)    ax[1].set_ylabel('Amplitude', fontsize=12)    ax[0].set_ylabel('Amplitude', fontsize=12)    ax[0].legend()    ax[1].legend()        plt.show()    # plt.subplots_adjust(hspace=0.5)#     if not figname:#         plt.savefig('fft_plus_power.png', dpi=300, bbox_inches='tight')#     else:#         plt.savefig(figname, dpi=300, bbox_inches='tight')#     plt.close('all')

plot_fft_plus_power(time, signal)# 对第二列的数值计算得到FFT能量谱序列

def plot_wavelet(time, signal, scales, waveletname='mexh', cmap=cmaps.GMT_panoply, title='Wavelet Transform (Power Spectrum) of signal', ylabel='Period (years)', xlabel='Time', figname=None):    # 丢进去一个时间参数x轴，序列值，伸缩尺度因子等    dt = time[1] - time[0]    #dt=1，对信号sidnal的第二列进行cwt操作，因为第一列是年份，进行变换的信号是一列，或者说一维    [coefficients, frequencies] = pywt.cwt(signal[:,0], scales, waveletname, dt)    # 进行连续小波小波变换，更多信息，包括dwt可以参看技术文档官网    power = (abs(coefficients)) ** 2    # 计算能量谱，为正数    period = 1. / frequencies
    scale0 = 8    numlevels = 10
    levels = [scale0]    for ll in range(1, numlevels):        # 1-9循环次数        scale0 *= 2        levels.append(scale0)
    contourlevels = np.log2(levels)    fig, ax = plt.subplots(figsize=(8, 6))    # 绘制三维图形，时间，周期和能量谱值    im = ax.contourf(time, np.log2(period), np.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**np.arange(np.ceil(np.log2(period.min())),                          np.ceil(np.log2(period.max())))    ax.set_yticks(np.log2(yticks))    ax.set_yticklabels(yticks)    ax.invert_yaxis()    # 倒转y轴    ylim = ax.get_ylim()


    
    # 自动获取y轴范围#     ax.set_ylim(ylim[0], -1)    cbar_ax = fig.add_axes([0.95, 0.15, 0.03, 0.7])    fig.colorbar(im, cax=cbar_ax, orientation="vertical")    plt.show()

plot_wavelet(time, signal, scales)