一维时间序列信号的奇异小波时频分析方法（Python）

最初的时频分析技术就是短时窗傅里叶变换STFT，由于时窗变短，可供分析的信号量减少，采用经典的谱估算方法引起的误差所占比重会增加。且该短时窗一旦选定．则在整个变换过程中其时窗长度是固定的。变换后的时频分辨率也即固定，这不利于低频和高频信号的检测。如果在时频分析中用长度变化的时窗来变换信号。就可更好地分析局部频率特征，这就是小波变换。根据测不准原理．小波变换的时间域分辨率和频率域分辨率不可能同时任意好，Stockwell提出的S变换是以Morlet小波为基本小波的连续小波变换的延伸．区别在于S变换中基本小波由简谐波与Gaussian函数乘积构成。基本小波中的简谐波在时间域仅作伸缩变换，Gaussian函数则进行伸缩和平移变换，而连续小波变换中简谐波与Gaussian函数同时进行伸缩和平移。此外，连续小波变换使用正交小波将信号从时间域分解到时间一尺度域，而S变换使用非正交小波将信号从时间域分解到时间一频率域。但由于S变换中的基本小波固定，无法根据地震信号本身选择基本小波。高静怀等Ⅲ提出了广义S变换．用带参数的调幅简谐波代替S变换中的基本小波，然后利用不同参数的基本小波作线性组合形成最终的基本小波。

有学者利用基于STFT时频分析技术对墨西哥湾密西西比河三角洲第四纪沉积环境进行了频谱成像。与传统的振幅属性相比，频谱成像后16Hz的峰值振幅属性图上河道的形态更清晰。与传统的相位属性相比，频谱成像后26Hz的相位属性图上河道的形态及断层均更清晰。

有学者给出了时频分析技术在储层含气性解释中的应用。对地震剖面进行频谱成像后，提取出了3个频率的振幅属性。从10Hz振幅属性剖面上可看到储层之下(约800ms处)有一个很强的反射能量．其也许为砂岩含气形成的“亮点”。但该“亮点”在20Hz振幅属性剖面上已经变弱，在30Hz振幅属性剖面上消失了．说明该“亮点”并不是由于砂岩含气形成的。

有学者提出了一种利用时频分析技术得到局部时频谱的谱反演方法，从而准确估算出随时间变化的地震子波，然后从地震数据中提取出时变子波并计算其反射系数，用于识别远小于地震分辨率的薄层。除了用于识别薄层外。该谱反演方法对刻画较小的上超或下超特征也很有效，这有助于分析准层序分布和沉积物搬运方向。

鉴于此，提出一种一维时间序列信号的奇异小波时频分析方法，时频分辨率得到了提高，程序采用Python编写。

import numpy as np
import waveletHelper # for wavelet functions

def morlet(Fc, Nc, Fs, norm):
    t = waveletHelper.getWaveletTimeRange(Fc, Nc, Fs) # get the time range where the wavelet is defined
    a = Fc*np.sqrt(2*np.pi)

    envelope = (a/Nc) * np.exp( - (t**2) * (a**2) * np.pi  / (Nc**2))
    wavelet = envelope * np.exp(1j * 2*np.pi*Fc * t)

    return waveletHelper.normalize(wavelet, envelope, norm, Fc)

def cwt(y, frange, Fs, baseCycle, norm, step):
    N, M = len(y), len(frange) # get N points from the input data and M frequency points for correct number of points operations    
    scalogram = np.zeros((M, N), dtype=float) # N points in time, M points in frequency

    for i in range(M):
        o = 1.0 + frange[i]/step # number of cycles increase factor per frequency
        w = morlet(frange[i], baseCycle*o, Fs, norm) # generate the wavelet
        lenW = len(w)
        scalogram[i,:] =  2 * np.abs( np.convolve(y, w, 'full')[lenW//2 -1 : -lenW//2] )**2 #apply convolution operation and correct for the resultant size
    
    return scalogram #rotate from time*frequency to frequency*time (helps heat
完整代码可通过知乎学术咨询获得：https://www.zhihu.com/consult/people/792359672131756032?isMe=1