1、從傅里葉變換到小波變換

傅里葉變換能夠?qū)⒁粋€(gè)信號(hào)從時(shí)域轉(zhuǎn)換為頻域，在轉(zhuǎn)換后的頻譜中，頻譜的峰值越大越尖，表示對(duì)應(yīng)頻率的信號(hào)就強(qiáng)度就越大。

傅里葉變換能夠處理不隨時(shí)間變化的平穩(wěn)信號(hào)，即它能告訴我們信號(hào)包含哪些頻段，但是不能告訴我們這個(gè)頻段是在信號(hào)的哪個(gè)時(shí)間段出現(xiàn)的。而生活中更多的是隨著時(shí)間變化而發(fā)生一定變化的非平穩(wěn)信號(hào)，而對(duì)于這些信號(hào)，傅里葉變換就難以進(jìn)行處理，因此演化出了可以處理非平穩(wěn)信號(hào)的小波變換。

如上圖所示，signal1和signal2兩個(gè)信號(hào)都是由4個(gè)不同頻率疊加得到的，只是疊加形式不一樣；對(duì)于signal1，四個(gè)不同頻率的信號(hào)直接進(jìn)行疊加，因?yàn)樾盘?hào)不隨時(shí)間變化而產(chǎn)生非平穩(wěn)變化，因此傅里葉變換可以輕松得到結(jié)果；但是對(duì)于信號(hào)signal2，四個(gè)不同頻率的信號(hào)以此在時(shí)間軸上展開，傅里葉變換就不能得出有效的結(jié)果了。

那么如何解決這個(gè)問題呢？短時(shí)傅里葉變換就出現(xiàn)了，想法很容易理解：如果將信號(hào)分割相同的n份，傅里葉變換找出特殊的信號(hào)在第幾個(gè)部分出現(xiàn)，那么我們就知道對(duì)應(yīng)的這個(gè)特殊信號(hào)頻段出現(xiàn)在原始信號(hào)的哪一段了。

但是這種方法也有問題，因?yàn)榘汛翱诘拇笮∽龅迷叫?，就越能知道信?hào)中某個(gè)頻率出現(xiàn)在哪里，但是自身的頻率值獲得的就越少。窗口的大小越大，我們對(duì)頻率值的了解就越多，而對(duì)時(shí)間的了解就越少。可以說是短時(shí)傅里葉變換不能夠兼顧頻域信息與時(shí)域信息。

因此短時(shí)傅里葉變換并不能很好解決這個(gè)問題，也可以說是引入了頻域信息與時(shí)域信息不兼容的問題。那么如何解決這個(gè)問題呢？這就引出了小波變換，它可以很好地處理這種動(dòng)態(tài)頻譜信息。小波變換不僅能夠告訴我們?cè)谛盘?hào)中有哪些頻率出現(xiàn)，而且能夠告訴我們?cè)谛盘?hào)的什么時(shí)候出現(xiàn)。小波變換的實(shí)現(xiàn)是通過不同的伸縮變換完成的。它能夠根據(jù)不同尺度的信號(hào)調(diào)節(jié)成不同大小的窗口，既能有效分析大尺度信號(hào)，也能有效分析小尺度信號(hào)，即既兼顧了時(shí)間信號(hào)，又兼顧了頻率信號(hào)。

下面用程序來說明這個(gè)問題：

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

t_n = 1
N = 200

xa = np.linspace(0, t_n, num=N)
xb = np.linspace(0, t_n / 4, num=N // 4)

def get_fft_values(y_values, N):
    xf = np.arange(len(y_values)) 
    f_values = xf[range(int(N / 2))]
    # f_values = np.linspace(0.0, 1.0,N // 2)
    fft_values_ = fft(y_values)
    fft_values = 2.0 / N * np.abs(fft_values_[0:N // 2])
    return f_values, fft_values


frequencies = [4,30, 60, 90]  
y1a, y1b = np.sin(2 * np.pi * frequencies[0] * xa), np.sin(2 * np.pi * frequencies[0] * xb)
y2a, y2b = np.sin(2 * np.pi * frequencies[1] * xa), np.sin(2 * np.pi * frequencies[1] * xb)
y3a, y3b = np.sin(2 * np.pi * frequencies[2] * xa), np.sin(2 * np.pi * frequencies[2] * xb)
y4a, y4b = np.sin(2 * np.pi * frequencies[3] * xa), np.sin(2 * np.pi * frequencies[3] * xb)

composite_signal1 = y1a + y2a + y3a + y4a  
composite_signal2 = np.concatenate([y1b, y2b, y3b, y4b])  

f_values1, fft_values1 = get_fft_values(composite_signal1, N)
f_values2, fft_values2 = get_fft_values(composite_signal2, N)

fig, axarr = plt.subplots(nrows=2, ncols=2, figsize=(8, 6))
axarr[0, 0].plot(xa, composite_signal1)
axarr[1, 0].plot(xa, composite_signal2)
axarr[0, 1].plot(f_values1, fft_values1)
axarr[1, 1].plot(f_values2, fft_values2)

plt.tight_layout()
plt.show()

赫茲 ，所有頻率在這段時(shí)間內(nèi)全部出現(xiàn)，右上角是這個(gè)信號(hào)的頻譜。
在下邊的兩個(gè)圖中，我們同樣能看到這4種頻率，但是第一種出現(xiàn)在1/4時(shí)間上，之后的頻率依次出現(xiàn)。并且我們能在右下面看到這個(gè)信號(hào)的頻譜圖。

值得注意的是兩個(gè)頻譜圖都包含了4個(gè)峰值，但是他們都不能明確表示出在信號(hào)的哪個(gè)時(shí)間段出現(xiàn)的。這就使得傅里葉變換不能區(qū)分這個(gè)兩個(gè)信號(hào)。

下面，我們用一張圖來說明小波變換與短時(shí)傅里葉變換的差異。
我們可以把豎狀長(zhǎng)方形當(dāng)作時(shí)間信息，把橫狀長(zhǎng)方形當(dāng)作頻域信息，短時(shí)傅里葉變換是通過統(tǒng)一大小窗口進(jìn)行變換，那么它在時(shí)間域和頻率域都只能有中等的分辨率。而右下角的小波變換卻可以根據(jù)信號(hào)需要變換長(zhǎng)方形的長(zhǎng)寬胖瘦，因而在時(shí)間域和頻率域上都能根據(jù)需要得到較好的分辨效果。

小波變換在時(shí)頻域變換特征為：
對(duì)于小頻率值，頻域分辨率高，時(shí)域分辨率低；
對(duì)于大頻率值，頻域分辨率低，時(shí)域分辨率高。

2、圖像化感受小波變換中的小波

傅里葉變換是通過一系列頻率不同的正弦波來擬合信號(hào)的。而小波變換是使用一系列稱為小波的函數(shù)(每個(gè)函數(shù)具有不同的尺度)擬合信號(hào)的。這里的小可以認(rèn)為是時(shí)間特征上的小。

由上圖我們可以看出，正弦波是在時(shí)間上無限延長(zhǎng)且保持周期變化，而小波則是在某一個(gè)時(shí)間點(diǎn)或時(shí)間段上有信號(hào)幅值，是局域波。因此小波不僅包含了頻率的信息，同時(shí)也包含了時(shí)間的信息。

小波是在時(shí)間域的函數(shù)，那么我們疊加不同時(shí)間段的小波，讓小波信號(hào)從研究信號(hào)的開始的時(shí)間點(diǎn)慢慢向結(jié)束時(shí)間點(diǎn)上移動(dòng)，即進(jìn)行卷積操作。在對(duì)原始小波操作完后，我們可以改變小波尺寸，然后繼續(xù)上面的操作。進(jìn)而就得到了研究信號(hào)的整個(gè)時(shí)域和頻域信息了。

3、小波變換族中的分類

傅里葉變換和小波變換的另一個(gè)不同的地方是小波中有很多不同的族。小波族之間是不同的，因?yàn)槊總€(gè)小波族在小波的緊湊性和平滑性方面都做出了不同的權(quán)衡。這意味著我們需要根據(jù)對(duì)應(yīng)的特征選擇一個(gè)特定的小波族。

通過pywt函數(shù)，可以看到小波變換中有多少個(gè)族類。
代碼及其結(jié)果顯示如下：

import pywt
print(pywt.families(short=False))

# 顯示的結(jié)果
['Haar', 'Daubechies', 'Symlets', 'Coiflets', 'Biorthogonal', 'Reverse biorthogonal', 
'Discrete Meyer (FIR Approximation)', 'Gaussian', 'Mexican hat wavelet', 'Morlet wavelet', 
'Complex Gaussian wavelets', 'Shannon wavelets', 'Frequency B-Spline wavelets', 'Complex Morlet wavelets']

小波之所以有很多的族類，對(duì)應(yīng)著不同的形狀、光滑度和緊密型，應(yīng)用于不同的場(chǎng)景和狀況，是因?yàn)樾〔ㄗ孱惖漠a(chǎn)生只需要滿足兩個(gè)數(shù)學(xué)條件：歸一化約束和正交化約束。即小波必須有：1、有限的能量；2、零均值。

有限能量意味著它在時(shí)間和頻率上是局域的，是可積的。零均值條件意味著小波在時(shí)域中均值為零，在時(shí)域中頻率為零時(shí)均值為零。這對(duì)于保證小波變換的可積性和計(jì)算小波變換的逆是必要的。

在此，我們可以查看不同小波族的圖像。下面圖像中第一行包含四個(gè)離散小波，第二行包含四個(gè)連續(xù)小波。

小波族圖像顯示代碼如下：

import pywt
from matplotlib import pyplot as plt

#print(pywt.families(short=False))
discrete_wavelets = ['db5', 'sym5', 'coif5', 'bior2.4']
continuous_wavelets = ['mexh', 'morl', 'cgau5', 'gaus5']

list_list_wavelets = [discrete_wavelets, continuous_wavelets]
list_funcs = [pywt.Wavelet, pywt.ContinuousWavelet]

fig, axarr = plt.subplots(nrows=2, ncols=4, figsize=(16,8))
for ii, list_wavelets in enumerate(list_list_wavelets):
    func = list_funcs[ii]
    row_no = ii
    for col_no, waveletname in enumerate(list_wavelets):
        wavelet = func(waveletname)
        family_name = wavelet.family_name
        biorthogonal = wavelet.biorthogonal
        orthogonal = wavelet.orthogonal
        symmetry = wavelet.symmetry
        if ii == 0:
            _ = wavelet.wavefun()
            wavelet_function = _[0]
            x_values = _[-1]
        else:
            wavelet_function, x_values = wavelet.wavefun()
        if col_no == 0 and ii == 0:
            axarr[row_no, col_no].set_ylabel("Discrete Wavelets", fontsize=16)
        if col_no == 0 and ii == 1:
            axarr[row_no, col_no].set_ylabel("Continuous Wavelets", fontsize=16)
        axarr[row_no, col_no].set_title("{}".format(family_name), fontsize=16)
        axarr[row_no, col_no].plot(x_values, wavelet_function)
        axarr[row_no, col_no].set_yticks([])
        axarr[row_no, col_no].set_yticklabels([])
plt.tight_layout()
plt.show()

4、小波族子類的分類

在每個(gè)小波族中，可以有許多不同的小波子類別屬于這個(gè)族?？梢酝ㄟ^系數(shù)的數(shù)量和分解的級(jí)別來區(qū)分小波的不同子類別。

下面我們以一個(gè)稱為“Daubechies”族為例來研究查看其子類。

上圖是多貝西小波族關(guān)于階數(shù)的分類情況，從第一列到第五列，分別對(duì)應(yīng)著db1到db5，當(dāng)然在此我們用的Daubechies小波可以高達(dá)20階，在此用5階來進(jìn)行研究和說明。

所用的階數(shù)表示消失力矩的個(gè)數(shù)，因此db3有3個(gè)消失力矩，db5有5個(gè)消失力矩。消失力矩的個(gè)數(shù)與小波的逼近階數(shù)和平滑度有關(guān)。如果一個(gè)小波有p個(gè)消失矩，那么它就可以近似為p - 1次多項(xiàng)式。

在選擇小波時(shí)，我們還可以指出分解的級(jí)別。默認(rèn)情況下，PyWavelets為輸入信號(hào)選擇可能的最大分解級(jí)別。分解的最大級(jí)別取決于輸入信號(hào)長(zhǎng)度和小波的長(zhǎng)度。(用pywt.dwt_max_level()來查看)

根據(jù)上圖可以看出，隨著消失力矩階數(shù)的增加，小波的多項(xiàng)式階數(shù)也增加，小波變得更加平滑。

對(duì)應(yīng)的代碼顯示如下：

import pywt
from matplotlib import pyplot as plt

db_wavelets = pywt.wavelist('db')[:5]
print(db_wavelets)
# ['db1', 'db2', 'db3', 'db4', 'db5']

fig, axarr = plt.subplots(ncols=5, nrows=5, figsize=(20,16))
fig.suptitle('Daubechies family of wavelets', fontsize=16)
for col_no, waveletname in enumerate(db_wavelets):
    wavelet = pywt.Wavelet(waveletname)
    no_moments = wavelet.vanishing_moments_psi
    family_name = wavelet.family_name
    for row_no, level in enumerate(range(1,6)):
        wavelet_function, scaling_function, x_values = wavelet.wavefun(level = level)
        axarr[row_no, col_no].set_title("{} - level {}\n{} vanishing moments\n{} samples".format(
            waveletname, level, no_moments, len(x_values)), loc='left')
        axarr[row_no, col_no].plot(x_values, wavelet_function, 'bD--')
        axarr[row_no, col_no].set_yticks([])
        axarr[row_no, col_no].set_yticklabels([])
plt.tight_layout()
plt.subplots_adjust(top=0.9)
plt.show()

5、連續(xù)小波變換與離散小波變換

小波變換有兩種不同的變換形式：連續(xù)小波變換和離散小波變換。

數(shù)學(xué)上，連續(xù)小波變換的表達(dá)式為:

其中Ψ(t)為連續(xù)母波，a是尺度因子，b為平移因子。尺度因子和平移因子的值是連續(xù)的，這意味著可能有無限多的小波。

離散小波變換和連續(xù)小波變換主要的區(qū)別是：離散小波變換對(duì)尺度和平移因子使用離散值?？s放因子a以2的冪增加(a=1,2,4…… )，平移因子b 整數(shù)值(b=1,2,3)增加。

6、離散小波變換–DWT是一組濾波器

在實(shí)際應(yīng)用中，離散小波變換總是被當(dāng)作濾波器使用。這意味著它被實(shí)現(xiàn)為高通濾波器和低通濾波器的級(jí)聯(lián)。這是因?yàn)闉V波器組是一種非常有效的方法，將一個(gè)信號(hào)分解成幾個(gè)子頻帶。

舉個(gè)例子，假設(shè)我們有一個(gè)頻率高達(dá)1000赫茲的信號(hào)。在第一個(gè)階段，我們把信號(hào)分成低頻部分和高頻部分，即0- 500hz和500- 1000hz。在第二階段，我們?nèi)〉皖l部分，再次將其分為兩部分:0-250赫茲和250-500赫茲。在第三階段，我們將0-250赫茲的部分分為0-125赫茲的部分和125-250赫茲的部分。這種情況會(huì)一直持續(xù)下去，直到我們達(dá)到了所需的精化水平，或者直到我們的樣本耗盡為止。

頻率劃分模式如下圖所示：

將一組信號(hào)均分為高頻信號(hào)和低頻信號(hào)后，如果進(jìn)一步再劃分，則將其低頻信號(hào)的頻率一分為二，再度產(chǎn)生一組新的高頻信號(hào)和低頻信號(hào)。上圖所示是五階小波劃分。

下面，我們用程序的形式來對(duì)一組波形進(jìn)行頻率的劃分。

import pywt
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 1, num=2048)
chirp_signal = np.sin(250 * np.pi * x ** 2)

fig, ax = plt.subplots(figsize=(6, 1))
ax.set_title("Original Chirp Signal: ")
ax.plot(chirp_signal)
plt.show()

data = chirp_signal
waveletname = 'sym5'

fig , axarr = plt.subplots(nrows=5, ncols=2, figsize=(6, 6))
for ii in range(5):
    (data, coeff_d) = pywt.dwt(data, waveletname)
    axarr[ii, 0].plot(data, 'r')
    axarr[ii, 1].plot(coeff_d, 'g')
    axarr[ii, 0].set_ylabel("Level {}".format(ii + 1), fontsize=14, rotation=90)
    axarr[ii, 0].set_yticklabels([])
    if ii == 0:
        axarr[ii, 0].set_title("Approximation coefficients", fontsize=14)
        axarr[ii, 1].set_title("Detail coefficients", fontsize=14)
    axarr[ii, 1].set_yticklabels([])
plt.tight_layout()
plt.show()

程序運(yùn)行后的結(jié)果顯示為：

原始chirp_signal信號(hào)顯示：

下圖是chirp_signal信號(hào)進(jìn)行sym5小波處理后(1 - 5級(jí))的近似系數(shù)和細(xì)節(jié)系數(shù)顯示，從1階到5階。左側(cè)一列是小波近似系數(shù)，右側(cè)是小波細(xì)節(jié)系數(shù)。

在小波中，離散小波變換(DWT)應(yīng)用函數(shù)pywt.dwt()；
DWT返回兩組系數(shù)：近似系數(shù)和細(xì)節(jié)系數(shù)；
近似系數(shù)表示小波變換的低通濾波器(平均濾波器)的輸出；
細(xì)節(jié)系數(shù)表示小波變換的高通濾波器(差分濾波器)的輸出；
通過對(duì)上一級(jí)小波變換的近似系數(shù)進(jìn)行小波變換，得到下一級(jí)小波變換；
每下一層，原始信號(hào)的采樣率也下降2倍。

7、連續(xù)小波變換實(shí)現(xiàn)狀態(tài)空間的可視化

為了更好理解小波變換尺度圖，讓我們將它與原始的時(shí)間序列數(shù)據(jù)及其傅里葉變換一起可視化。

下面程序所用的數(shù)據(jù)來自一段ppg信號(hào)，讀者可以自己選取一段任意的信號(hào)作為dataset輸入即可。

import matplotlib.pyplot as plt
import numpy as np
import pywt
from scipy.fftpack import fft
import os
from scipy import signal


def plot_wavelet(time, signal, scales,
                 waveletname='cmor',
                 cmap=plt.cm.seismic,
                 title='Wavelet Transform (Power Spectrum) of signal',
                 ylabel='Period (years)',
                 xlabel='Time'):
    dt = time[1] - time[0]
    [coefficients, frequencies] = pywt.cwt(signal, scales, waveletname, dt)
    power = (abs(coefficients)) ** 2
    period = 1. / frequencies
    levels = [0.0625, 0.125, 0.25, 0.5, 1, 2, 4, 8]
    contourlevels = np.log2(levels)

    fig, ax = plt.subplots(figsize=(15, 10))
    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()
    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")
    # plt.show()

def get_ave_values(time, signal, average_over):
    # 進(jìn)行均值濾波去噪
    smooth_list=[]
    myfilter=[]
    filter_window=average_over
    for i in range(len(signal)):
        if len(myfilter) < filter_window:  # 與濾波窗口大小作比較
            myfilter.append( signal[i])
        else:
            smooth_data = np.mean(myfilter)
            myfilter.pop(0)  # 減去第一個(gè)數(shù)據(jù)，隨后加一個(gè)新的數(shù)，形成一個(gè)動(dòng)態(tài)的均值濾波器
            myfilter.append(signal[i])
            # raw_list.append(raw_data)
            smooth_list.append(smooth_data)
    return time,smooth_list

def plot_signal_plus_average(time, signal, average_over=5):
    fig, ax = plt.subplots(figsize=(15, 3))
    time_ave, signal_ave = get_ave_values(time, signal, average_over)
    ax.plot(time, signal, label='signal')
    ax.plot(signal_ave, label='time average (n={})'.format(5))
    ax.set_xlim([time[0], time[-1]])
    ax.set_ylabel('Signal Amplitude', fontsize=18)
    ax.set_title('Signal + Time Average', fontsize=18)
    ax.set_xlabel('Time', fontsize=18)
    ax.legend()
    # plt.show()


def get_fft_values(y_values,N):
    xf = np.arange(len(y_values))  # 頻率
    f_values = xf[range(int(N / 2))]
    fft_values_ = fft(y_values)
    fft_values = 2.0 / N * np.abs(fft_values_[0:N // 2])
    return f_values, fft_values


def plot_fft_plus_power(time, signal):
    N = len(signal)
    fig, ax = plt.subplots(figsize=(15, 3))
    variance = np.std(signal) ** 2
    time_ave, signal_ave = get_ave_values(time, signal, average_over=5)
    f_values, fft_values = get_fft_values(signal_ave, N)
    fft_power = variance * abs(fft_values) ** 2  # FFT power spectrum
    ax.plot(f_values, fft_values, 'r-', label='Fourier Transform')
    ax.plot(f_values, fft_power, 'k--', linewidth=1, label='FFT Power Spectrum')
    ax.set_xlabel('Frequency [Hz / year]', fontsize=18)
    ax.set_ylabel('Amplitude', fontsize=18)
    ax.legend()
    # plt.show()

fname = os.path.join("D:\\a_user_file\\data\\1_Pulse Transit Time PPG Dataset", "s1_run.csv")
dataset = []

with open(fname, "r") as f1:
    data = f1.read()
    lines = data.split("\n")
    #lines = lines[1:len(lines) - 1]
    lines = lines[1:3000]
    num = []
    for i in range(len(lines)):
        line_i = lines[i].split(",")
        num.append(int(line_i[4]))
    print(len(num))  # 245274
    for i in range(len(num)):
        if i%5==0:
            dataset.append(int(num[i]))
    print(len(dataset))

time = np.arange(len(dataset))
signal = dataset
scales = np.arange(1, 128)

plot_signal_plus_average(time, signal)
plot_fft_plus_power(time, signal)
plot_wavelet(time, signal, scales)
plt.show()

最后顯示的結(jié)果為：

8、小波分解重構(gòu)信號(hào)

小波是如何將一個(gè)信號(hào)分解成它的子頻帶，并重新構(gòu)造原始信號(hào)呢？這一節(jié)將詳細(xì)介紹：

PyWavelets 提供了兩種解析信號(hào)的方法：
(1)用pywt.dwt()函數(shù)來檢索近似系數(shù)。然后對(duì)檢索到的系數(shù)應(yīng)用DWT以獲得第二級(jí)系數(shù)，并繼續(xù)此過程，直到達(dá)到所需的分解級(jí)別；
(2)用pywt.wavedec()函數(shù)直接和檢索所有的細(xì)節(jié)系數(shù)達(dá)到某種程度上n。這個(gè)函數(shù)作為輸入原始信號(hào)和n水平并返回一組近似系數(shù)(第n個(gè)水平)和n組細(xì)節(jié)系數(shù)(1到n級(jí))。

8.1、pywt.dwt()函數(shù)解構(gòu)重構(gòu)信號(hào)

import pywt
import numpy as np
import matplotlib.pyplot as plt
import os

# fname是存放數(shù)據(jù)的目錄
fname = os.path.join("D:\\a_user_file\\data\\1_Pulse Transit Time PPG Dataset", "s1_run.csv")

with open(fname, "r") as f1:
    data = f1.read()
    lines = data.split("\n")
    #lines = lines[1:len(lines) - 1]
    lines = lines[1:3000]
    data_set = []
    for i in range(len(lines)):
        line_i = lines[i].split(",")
        data_set.append(int(line_i[4]))
    # print(len(data_set))  
    
signal=data_set[:3000]
(cA1, cD1) = pywt.dwt(signal, 'db2', 'smooth')
reconstructed_signal = pywt.idwt(cA1, cD1, 'db2', 'smooth')

plt.subplot(211)
plt.plot(signal, label='signal')
plt.xlabel('x-采樣點(diǎn)', fontproperties='SimHei')
plt.ylabel('y-幅值', fontproperties='SimHei')
plt.title('原始波形', fontproperties='STLITI', fontsize=10)
plt.subplot(212)
plt.plot(reconstructed_signal, label='reconstructed signal', linestyle='--')
plt.xlabel('x-采樣點(diǎn)', fontproperties='SimHei')
plt.ylabel('y-幅值', fontproperties='SimHei')
plt.title('合成波形', fontproperties='STLITI', fontsize=10)
plt.show()

運(yùn)行后的結(jié)果顯示圖為：

上面圖像中實(shí)線顯示為原始數(shù)據(jù)圖像，虛線為pywt.dwt()函數(shù)解析后構(gòu)造出來的信號(hào)圖像顯示，從幅值和輪廓上都可以看出pywt.dwt()函數(shù)可以比較完美得構(gòu)造出原始信號(hào)。其中pywt.dwt()函數(shù)更多用于低級(jí)別系數(shù)的重構(gòu)。

8.2、pywt.wavedec()函數(shù)解構(gòu)重構(gòu)信號(hào)

用9.1所用的數(shù)據(jù)集進(jìn)行檢驗(yàn)并對(duì)比，將中間的信號(hào)解構(gòu)重構(gòu)函數(shù)修改后，程序如下所示：

import pywt
import numpy as np
import matplotlib.pyplot as plt
import os

# fname是存放數(shù)據(jù)的目錄
fname = os.path.join("D:\\a_user_file\\data\\1_Pulse Transit Time PPG Dataset", "s1_run.csv")

with open(fname, "r") as f1:
    data = f1.read()
    lines = data.split("\n")
    #lines = lines[1:len(lines) - 1]
    lines = lines[1:3000]
    data_set = []
    for i in range(len(lines)):
        line_i = lines[i].split(",")
        data_set.append(int(line_i[4]))
    # print(len(data_set))  
    
signal=data_set[:3000]
coeffs = pywt.wavedec(signal, 'db2', level=8)
reconstructed_signal = pywt.waverec(coeffs, 'db2')

plt.subplot(211)
plt.plot(signal, label='signal')
plt.xlabel('x-采樣點(diǎn)', fontproperties='SimHei')
plt.ylabel('y-幅值', fontproperties='SimHei')
plt.title('原始波形', fontproperties='STLITI', fontsize=10)
plt.subplot(212)
plt.plot(reconstructed_signal, label='reconstructed signal', linestyle='--')
plt.xlabel('x-采樣點(diǎn)', fontproperties='SimHei')
plt.ylabel('y-幅值', fontproperties='SimHei')
plt.title('合成波形', fontproperties='STLITI', fontsize=10)
plt.show()

運(yùn)行后顯示的圖像為：

可見，pywt.wavedec()函數(shù)也能夠很好得將信號(hào)重構(gòu)出來，其中pywt.wavedec()函數(shù)更多用于高級(jí)別系數(shù)的重構(gòu)。

9、使用離散小波變換去除噪聲信號(hào)

根據(jù)上面內(nèi)容，我們知道可以將信號(hào)分解為近似(低通)和細(xì)節(jié)(高通)系數(shù)。并且可以用這些系數(shù)重建信號(hào)，我們將得到原始信號(hào)。

如果我們?cè)谥亟ǖ倪^程中丟掉一些細(xì)節(jié)系數(shù)會(huì)怎么樣呢？因?yàn)榧?xì)節(jié)系數(shù)代表的是信號(hào)中的高頻部分，我們只需要過濾掉這部分頻譜，就相當(dāng)于可以過濾掉高頻的噪聲信號(hào)。

此處我們可以使用pywt.threshold()來刪除細(xì)節(jié)系數(shù)，它可以刪除了高于給定閾值的系數(shù)值。在進(jìn)行解構(gòu)信號(hào)后，把一些系數(shù)設(shè)為零，然后重新構(gòu)造，我們可以從信號(hào)中去除高頻噪聲。

此處，所用的pywt.threshold()函數(shù)來刪除細(xì)節(jié)系數(shù)的代碼顯示如下：

def lowpassfilter(signal, thresh = 0.63, wavelet="db4"):
    thresh = thresh*np.nanmax(signal)
    coeff = pywt.wavedec(signal, wavelet, mode="per" )
    coeff[1:] = (pywt.threshold(i, value=thresh, mode="soft" ) for i in coeff[1:])
    reconstructed_signal = pywt.waverec(coeff, wavelet, mode="per" )
    return reconstructed_signal

rec = lowpassfilter(signal, 0.4)

10、使用離散小波變換進(jìn)行信號(hào)分類

之前的內(nèi)容可以看到如何使用連續(xù)小波變換和卷積神經(jīng)網(wǎng)絡(luò)來對(duì)信號(hào)進(jìn)行分類的。其實(shí)用離散小波變換也可以進(jìn)行信號(hào)的分類的。

離散小波變換分類思想如下：離散小波變換用于將一個(gè)信號(hào)分割成不同的子頻帶，盡可能多的子頻帶或盡可能多的子頻帶。如果不同類型的信號(hào)表現(xiàn)出不同的頻率特性，這種行為上的差異必須表現(xiàn)在一個(gè)頻率子帶中。因此，如果我們從每個(gè)子帶生成特征，并將這些特征集合作為分類器的輸入(Random Forest, Gradient boost, Logistic Regression等)，并利用這些特征訓(xùn)練分類器，那么分類器應(yīng)該能夠區(qū)分不同類型的信號(hào)。即：將信號(hào)進(jìn)行離散小波變換時(shí)，可以分解它的頻率子帶。并從每個(gè)子帶中提取特征，作為分類器的輸入。

那么，可以從每個(gè)子帶的數(shù)值集中生成什么樣的特性呢?當(dāng)然，這在很大程度上取決于信號(hào)的類型和應(yīng)用程序。但是一般來說，依然有一些最常用的特性。

(1)自回歸模型系數(shù)值;
(2)熵值:熵值可以用來衡量信號(hào)的復(fù)雜度;
(3)統(tǒng)計(jì)特征如：方差、標(biāo)準(zhǔn)差、均值、中值、第25百分位值、第75個(gè)百分位值、均方根值、導(dǎo)數(shù)的均值、過零率：即信號(hào)通過y = 0的次數(shù)。

def calculate_entropy(list_values):
    counter_values = Counter(list_values).most_common()
    probabilities = [elem[1]/len(list_values) for elem in counter_values]
    entropy=scipy.stats.entropy(probabilities)
    return entropy

def calculate_statistics(list_values):
    n5 = np.nanpercentile(list_values, 5)
    n25 = np.nanpercentile(list_values, 25)
    n75 = np.nanpercentile(list_values, 75)
    n95 = np.nanpercentile(list_values, 95)
    median = np.nanpercentile(list_values, 50)
    mean = np.nanmean(list_values)
    std = np.nanstd(list_values)
    var = np.nanvar(list_values)
    rms = np.nanmean(np.sqrt(list_values**2))
    return [n5, n25, n75, n95, median, mean, std, var, rms]

def calculate_crossings(list_values):
    zero_crossing_indices = np.nonzero(np.diff(np.array(list_values) > 0))[0]
    no_zero_crossings = len(zero_crossing_indices)
    mean_crossing_indices = np.nonzero(np.diff(np.array(list_values) > np.nanmean(list_values)))[0]
    no_mean_crossings = len(mean_crossing_indices)
    return [no_zero_crossings, no_mean_crossings]

def get_features(list_values):
    entropy = calculate_entropy(list_values)
    crossings = calculate_crossings(list_values)
    statistics = calculate_statistics(list_values)
    return [entropy] + crossings + statistics

上面代碼可以看到有以下四個(gè)函數(shù)：(1)、計(jì)算輸入信號(hào)熵值的函數(shù)；(2)、用來計(jì)算一些統(tǒng)計(jì)數(shù)據(jù)的函數(shù)，如幾個(gè)百分位數(shù)、平均值、標(biāo)準(zhǔn)差、方差等；(3)、計(jì)算過零率和平均過零率的函數(shù)；(4)將這三個(gè)函數(shù)的結(jié)果結(jié)合起來的函數(shù)。

最后一個(gè)get_features()函數(shù)為任何值列表返回一組12個(gè)特性。因此，如果將一個(gè)信號(hào)分解成10個(gè)不同的子帶，我們?yōu)槊總€(gè)子帶生成特征，我們將得到每個(gè)信號(hào)的10*12 = 120個(gè)特征。

參考文本：http://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/文章來源地址http://www.zghlxwxcb.cn/news/detail-813257.html

到了這里，關(guān)于從傅里葉變換到小波變換詳細(xì)解釋(含代碼)的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！