傅里葉變換


傅里葉變換是用三角函數表示目標函數,傅里葉變換廣泛的應用在信號處理、偏微分方程、熱力學、概率統計等領域:大到天體觀測,小到我們手機中圖片、音頻應用等,沒有傅里葉變換就沒有如今豐富多彩的信息化時代。在人工智能領域中,可利用傅里葉變換證明中心極限定理,而中心極限定理是概率學最重要的基石;傅里葉變換本質是將時域的信息匯總到頻域中,當兩組數據的傅里葉變換結果相同時,稱為兩者依概率收斂。本章介紹傅里葉變換推導過程並結合代碼實現傅里葉變換,按數學推導的離散傅里葉變換算法復雜度較高,本章最后介紹快速傅里葉變換(FFT)的實現原理。

一、傅里葉變換原理

1.1 目標函數f(x)為周期2π的函數

前面的章節中曾用多項式擬合目標函數,傅里葉變換是利用三角函數擬合函數,正弦、余弦函數有以下性質:

正交.png (1)

三角函數集合正交基礎.png在[-π,π]可以看成函數空間一組正交基,目標函數f(x)可寫成傅里葉級數:

三角函數形式1.png (1.1)

an、bn都是待定系數,也表示f(x)在正交基上的坐標,(1.1)式兩邊乘以cos mx ,m為一個特定系數坐標,求積分有:

推導1.png

由此可得an:

an.png    (1.3)

同樣可得bn:

bn.png      (1.4)

上面函數f(x)周期為2π,如果f(x)是周期為2T,可用線性仿射將f(x)的定義域x投射到t上,t的周期是2π:

放射.png

此時有變換1.png,變元.png,從上圖可見,Φ(t)是在[-∞,∞]一個周期為2π函數,有:

pt.png

變換1.png轉化為4png.png代入上式可得:

周期t.png (1.5)

再利用積分求系數an、bn:

ttt.png  (1.6)

1.2 f(x)非周期函數

f(x)不是周期函數時,可理解為f(x)的周期為+∞,不妨先考慮周期為2T函數fT(x),T為無限大取極限后,可將fT(x)延拓到周期為+∞函數f(x),先利用歐拉公式將公式(1.5)變為復數形式:

歐拉公式.png

頻率.png公式(1.5)可表示成:

傅里葉級數轉化為復數形式.png     (1.7)

an+ibn和an-ibn是一對共軛復數,幅值相同,相位角在x軸上對稱,設cn=an-ibn,利用公式(1.6)得:

cn.png    (1.7.1)

cn的共軛共軛2.png且有:

c2n.png  (1.7.2)

公式1.7中n的取值范圍是從1到無窮,引入cn后n的取值范圍是(-∞,+∞),1.7式寫成下面的形式:

復數技術形式.png   (1.7.3)

1.7.1代入上式后得到fT(x):

FT21.png   (1.8)

接下來利用式(1.8)將fT(x)周期延拓至∞可得目標函數f(x):

TW.png    (1.8.1)

式(1.8.1)代入(1.8)得到f(x)極限形式的函數:

fft21.png

當T為+∞時,wn.png,f(x)此時可表示成積分形式:

傅里葉變換.png (1.9)

上式中傅里葉變換積分形式.png稱為原函數f(x)的傅里葉變換,計算積分過程中ω視為符號常量,傅里葉變換積分形式.png積分結果是含自變量ω的函數,傅里葉變換可用下面的表示方式:

傅里葉變換形式1.png (1.10)

而F(f)再經過一次積分后可還原為f(x),這稱為傅里葉逆變換:

傅里葉逆變換.png

f(x)同樣稱為時域函數,比如記錄每天雨量的變化、每小時道路車輛的流量、汽車每年的里程數等,時域函數是日常生活中大量接觸到的函數模型,而1.10式將f(x)經過傅里葉變換后得到是頻域的信息,頻域單位是赫茲,反應出原函數周期性方面的信息,時域與頻域觀察世界視角不同,可以用下圖來加強對兩者的認識。

時域與頻域.png

二、程序實現傅里葉變換

2.1 離散傅里葉變換

代碼中利用式(1.10)的離散形式實現傅里葉變換,設在定義域內選擇N個數據實現f(x)函數離散化,f(x)可表示成:

傅里葉變換離散形式.png (2)

如果f(x)表示的是一個信號,所選擇的時間段為一秒即單位時間,那么N稱為采樣頻率,由采樣定理可知,采樣頻率至少是信號最大頻率的2倍,就傅里葉變換而言,采樣頻率最好是2的整數倍,有了這些設定后公式(2)可表示為:

離散共1.png    (2.1)

上式用12.png將單位時間一秒分成N份,在傅里葉變換中也稱歸一化,2.1式中ω代表角速度,通常將角速度處理成頻率形式,頻率h與角速度的關系有:

匹配.png

上式代入2.1后有:

xh.png  (2.2)

X(h)是單位為赫茲的頻域函數,再回頭看公式1.7.3,參數cn中的n取值范圍是(-∞,+∞),也就是說原公式中是允許有負角速度或者說負頻率存在的,而(2.2)式中n是非負整數,代表着客觀物理世界中沒有所謂的負頻率,這是數學與物理的差異性導致的,(2.2)引入負頻率才能滿足傅里葉變換公式,幸運的是正頻率與負頻率是共軛的關系,即cn與cn好.png關系,兩者的幅值是一樣的,傅里葉離散變換的目標是得到每個頻率對應幅值,所以將(2.2)乘以2即可得到與1.7.3式一樣的定義。

頻譜2.png      (2.3)

下面代碼利用公式2.3實現傅里葉離散變換,目標函數(原始信號)為:

testf.png

原始信號含有頻率為180、390、600Hz的信號,采用頻率為2048Hz,對應幅值分別為7、1.5、5.1。

import numpy as np
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from matplotlib.pylab import mpl

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 顯示中文
mpl.rcParams['axes.unicode_minus'] = False  # 顯示負號
Samplingfq=2048
def loadsignal(N  ):
    # 采樣點選擇1400個,因為設置的信號頻率分量最高為600Hz,根據采樣定理知采樣頻率要大於信號頻率2倍,所以這里設置采樣頻率為1400Hz(即一秒內有1400個采樣點)
    x = np.linspace(0, 1, N)
    # 設置需要采樣的信號,頻率分量有180,390和600
    y = 7 * np.sin(2 * np.pi * 180 * x) + 1.5 * np.sin(2 * np.pi * 390 * x) + 5.1 * np.sin(2 * np.pi * 600 * x)
    return  x,y
 

def FourierTransform():
    x, y = loadsignal(Samplingfq)
    ffts=[]
    for i in range(Samplingfq):
        real=0
        imag=0
        for t in range(Samplingfq):
            real=real+y[t]* np.cos(2*np.pi*i*x[t])*2/Samplingfq
            imag=imag+y[t]*  np.sin(2*np.pi*i*x[t])*-2/Samplingfq
        ffts.append([i,np.sqrt(real**2+imag**2)  ])
    ffts=np.array(ffts)
    plt.subplot(211)
    plt.plot(range(int(Samplingfq/2)),  ffts[:int(Samplingfq/2),1]  , 'b')
    plt.title('歸一化單邊', fontsize=10, color='#F08080')
    plt.subplot(212)
    plt.plot(range(Samplingfq),  ffts[:, 1]  , 'b')
    plt.title('歸一化雙邊', fontsize=10, color='#F08080')
    plt.show()


if __name__ == '__main__':
    FourierTransform()
 for i in range(Samplingfq):
        real=0
        imag=0
        for t in range(Samplingfq):
            real=real+y[t]* np.cos(2*np.pi*i*x[t])*2/Samplingfq
            imag=imag+y[t]*  np.sin(2*np.pi*i*x[t])*-2/Samplingfq

上面代碼片段即為對公式2.3應用,傅里葉離散變換后得到頻率圖如下圖所示:

Figure_1.png

2.2 快速傅里葉變換FFT

上面的示例程序復雜度是n2,n對應於采樣頻率,當目標函數或原始信號的最大頻率較高時,根據采樣定理,采樣頻率至少為最大頻率的2倍,這就造成n2值非常大,算法復雜度較高影響傅里葉變換的使用,利用快速傅里葉變換(FFT)可將復雜度降至復雜度.png,隨着采樣頻率的升高,FFT效率提升越來越明顯,FFT是信號傳輸、圖像分析、頻譜分析、數據壓縮最重要的數據算法之一,FFT出現給傅里葉變換帶來了旺盛的生命力。

單位時間內采樣頻率為N,0到N-1之間任意時刻t信號值用x(t)表示,2.3式用下式表示:

信號表示1.png (2.4)

上式中k代表0到N-1之間任意一個頻率,且有0<=k<N-1,引入一個新的符號變量:

單位圓的根.png

wnk.png是單位圓xN=1的一個根,wnk.png有以下性質:

wnk性質.png

利用歐拉公式、三角函數周期性即可證明這兩個性質,引入wnk.png后2.4式可寫為:

信號形式1.png    (2.5)

X(k)內含一個N-1次多項式,系數是x(0)、x(1)、x(2)....x(N-1) ,針對該多項式有以下性質:

多項式.png

在引入另外兩個多項式A1(x)、A2(x):

多項式2.png

三者之間關系為:

多項式44.png    (2.6)

wnk.png代入到2.6式中,並結合wnk.png自身的性質,當k<N/2時有:

GONGSI1.png   (2.6.1)

k+N/2代入2.6.1式有:

G21.png   (2.6.2)

由式2.6.1和式2.6.2可知,得到NK.png后,即可知道WKNN12.png值,而WKNN12.png再經過簡單處理后即可得到傅里葉變換后的值:

31.png     (2.7)

下面以一個簡單信號為例,在單位時間內采樣頻率為8,8個采樣點信號分別為0,1,2,3,4,5,6,7。對這組信號采用快速傅里葉變換有以下過程。簡單起見,下圖中數值即代表信號編號也代表信號的數值大小:

1624863153592068447.png

以求a81.png為例,root1.png,a81.png可分解為:

第一層.png

這個分解對應於樹狀圖的第一層,其中a141.png對應系數是0、2、4、6,A241.png對應的系數是1、3、5、7;a141.png可以理解為一個新的傅里葉變換:

第二次.png

余下文章請轉至鏈接:傅里葉變換

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM