|
目錄 一、一些關(guān)鍵概念的引入 1、離散傅里葉變換（DFT） 2、快速傅里葉變換（FFT） 3、采樣頻率以及采樣定理 4、如何理解采樣定理？ 二、使用scipy包實(shí)現(xiàn)快速傅里葉變換 1、產(chǎn)生原始信號(hào)——原始信號(hào)是三個(gè)正弦波的疊加 2、快速傅里葉變換 3、FFT的原始頻譜 4、將振幅譜進(jìn)行歸一化和取半處理 三、完整代碼 一、一些關(guān)鍵概念的引入 1、離散傅里葉變換（DFT） 離散傅里葉變換(discrete Fourier transform) 傅里葉分析方法是信號(hào)分析的最基本方法，傅里葉變換是傅里葉分析的核心，通過它把信號(hào)從時(shí)間域變換到頻率域，進(jìn)而研究信號(hào)的頻譜結(jié)構(gòu)和變化規(guī)律。但是它的致命缺點(diǎn)是：計(jì)算量太大，時(shí)間復(fù)雜度太高，當(dāng)采樣點(diǎn)數(shù)太高的時(shí)候，計(jì)算緩慢，由此出現(xiàn)了DFT的快速實(shí)現(xiàn)，即下面的快速傅里葉變換FFT。 2、快速傅里葉變換（FFT） 計(jì)算量更小的離散傅里葉的一種實(shí)現(xiàn)方法。詳細(xì)細(xì)節(jié)這里不做描述。 3、采樣頻率以及采樣定理 采樣頻率：采樣頻率，也稱為采樣速度或者采樣率，定義了每秒從連續(xù)信號(hào)中提取并組成離散信號(hào)的采樣個(gè)數(shù)，它用赫茲（Hz）來表示。采樣頻率的倒數(shù)是采樣周期或者叫作采樣時(shí)間，它是采樣之間的時(shí)間間隔。通俗的講采樣頻率是指計(jì)算機(jī)每秒鐘采集多少個(gè)信號(hào)樣本。 采樣定理：所謂采樣定理 ，又稱香農(nóng)采樣定理，奈奎斯特采樣定理，是信息論，特別是通訊與信號(hào)處理學(xué)科中的一個(gè)重要基本結(jié)論。采樣定理指出，如果信號(hào)是帶限的，并且采樣頻率高于信號(hào)帶寬的兩倍，那么，原來的連續(xù)信號(hào)可以從采樣樣本中完全重建出來。 定理的具體表述為：在進(jìn)行模擬/數(shù)字信號(hào)的轉(zhuǎn)換過程中，當(dāng)采樣頻率fs大于信號(hào)中最高頻率fmax的2倍時(shí),即 fs>2*fmax 采樣之后的數(shù)字信號(hào)完整地保留了原始信號(hào)中的信息，一般實(shí)際應(yīng)用中保證采樣頻率為信號(hào)最高頻率的2.56～4倍；采樣定理又稱奈奎斯特定理。 4、如何理解采樣定理？ 在對(duì)連續(xù)信號(hào)進(jìn)行離散化的過程中，難免會(huì)損失很多信息，就拿一個(gè)簡(jiǎn)單地正弦波而言，如果我1秒內(nèi)就選擇一個(gè)點(diǎn)，很顯然，損失的信號(hào)太多了，光著一個(gè)點(diǎn)我根本不知道這個(gè)正弦信號(hào)到底是什么樣子的，自然也沒有辦法根據(jù)這一個(gè)采樣點(diǎn)進(jìn)行正弦波的還原，很明顯，我采樣的點(diǎn)越密集，那越接近原來的正弦波原始的樣子，自然損失的信息越少，越方便還原正弦波。故而 采樣定理說明采樣頻率與信號(hào)頻率之間的關(guān)系，是連續(xù)信號(hào)離散化的基本依據(jù)。 它為采樣率建立了一個(gè)足夠的條件，該采樣率允許離散采樣序列從有限帶寬的連續(xù)時(shí)間信號(hào)中捕獲所有信息。 二、使用scipy包實(shí)現(xiàn)快速傅里葉變換 本節(jié)不會(huì)說明FFT的底層實(shí)現(xiàn)，只介紹scipy中fft的函數(shù)接口以及使用的一些細(xì)節(jié)。 1、產(chǎn)生原始信號(hào)——原始信號(hào)是三個(gè)正弦波的疊加 import numpy as npfrom scipy.fftpack import fft,ifftimport matplotlib.pyplot as pltfrom matplotlib.pylab import mpl mpl.rcParams['font.sans-serif'] = ['SimHei'] #顯示中文mpl.rcParams['axes.unicode_minus']=False #顯示負(fù)號(hào) #采樣點(diǎn)選擇1400個(gè)，因?yàn)樵O(shè)置的信號(hào)頻率分量最高為600赫茲，根據(jù)采樣定理知采樣頻率要大于信號(hào)頻率2倍，所以這里設(shè)置采樣頻率為1400赫茲（即一秒內(nèi)有1400個(gè)采樣點(diǎn)，一樣意思的）x=np.linspace(0,1,1400) #設(shè)置需要采樣的信號(hào)，頻率分量有200，400和600y=7*np.sin(2*np.pi*200*x) + 5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x) 這里原始信號(hào)的三個(gè)正弦波的頻率分別為，200Hz、400Hz、600Hz,最大頻率為600赫茲。根據(jù)采樣定理，fs至少是600赫茲的2倍，這里選擇1400赫茲，即在一秒內(nèi)選擇1400個(gè)點(diǎn)。 原始的函數(shù)圖像如下： plt.figure()plt.plot(x,y) plt.title('原始波形') plt.figure()plt.plot(x[0:50],y[0:50]) plt.title('原始部分波形（前50組樣本）')plt.show()由圖可見，由于采樣點(diǎn)太過密集，看起來不好看，我們只顯示前面的50組數(shù)據(jù)，如下： 2、快速傅里葉變換 其實(shí)scipy和numpy一樣，實(shí)現(xiàn)FFT非常簡(jiǎn)單，僅僅是一句話而已，函數(shù)接口如下： from scipy.fftpack import fft,ifft from numpy import fft,ifft 其中fft表示快速傅里葉變換，ifft表示其逆變換。具體實(shí)現(xiàn)如下： fft_y=fft(y) #快速傅里葉變換print(len(fft_y))print(fft_y[0:5])'''運(yùn)行結(jié)果如下：1400[-4.18864943e-12+0.j 9.66210986e-05-0.04305756j 3.86508070e-04-0.08611996j 8.69732036e-04-0.12919206j 1.54641157e-03-0.17227871j]''' 我們發(fā)現(xiàn)以下幾個(gè)特點(diǎn)： （1）變換之后的結(jié)果數(shù)據(jù)長(zhǎng)度和原始采樣信號(hào)是一樣的 （2）每一個(gè)變換之后的值是一個(gè)復(fù)數(shù)，為a+bj的形式，那這個(gè)復(fù)數(shù)是什么意思呢？ 我們知道，復(fù)數(shù)a+bj在坐標(biāo)系中表示為（a,b），故而復(fù)數(shù)具有模和角度，我們都知道快速傅里葉變換具有 “振幅譜”“相位譜”，它其實(shí)就是通過對(duì)快速傅里葉變換得到的復(fù)數(shù)結(jié)果進(jìn)一步求出來的， 那這個(gè)直接變換后的結(jié)果是不是就是我需要的，當(dāng)然是需要的，在FFT中，得到的結(jié)果是復(fù)數(shù)， （3）FFT得到的復(fù)數(shù)的模（即絕對(duì)值）就是對(duì)應(yīng)的“振幅譜”，復(fù)數(shù)所對(duì)應(yīng)的角度，就是所對(duì)應(yīng)的“相位譜”，現(xiàn)在可以畫圖了。 3、FFT的原始頻譜 N=1400x = np.arange(N) # 頻率個(gè)數(shù) abs_y=np.abs(fft_y) # 取復(fù)數(shù)的絕對(duì)值，即復(fù)數(shù)的模(雙邊頻譜)angle_y=np.angle(fft_y) #取復(fù)數(shù)的角度 plt.figure()plt.plot(x,abs_y) plt.title('雙邊振幅譜（未歸一化）') plt.figure()plt.plot(x,angle_y) plt.title('雙邊相位譜（未歸一化）')plt.show()顯示結(jié)果如下： 注意：我們?cè)诖颂巸H僅考慮“振幅譜”，不再考慮相位譜。 我們發(fā)現(xiàn)，振幅譜的縱坐標(biāo)很大，而且具有對(duì)稱性，這是怎么一回事呢？ 關(guān)鍵：關(guān)于振幅值很大的解釋以及解決辦法——?dú)w一化和取一半處理 比如有一個(gè)信號(hào)如下： Y=A1+A2*cos(2πω2+φ2）+A3*cos(2πω3+φ3）+A4*cos(2πω4+φ4） 經(jīng)過FFT之后，得到的“振幅圖”中， 第一個(gè)峰值（頻率位置）的模是A1的N倍，N為采樣點(diǎn)，本例中為N=1400，此例中沒有，因?yàn)樾盘?hào)沒有常數(shù)項(xiàng)A1 第二個(gè)峰值（頻率位置）的模是A2的N/2倍，N為采樣點(diǎn)， 第三個(gè)峰值（頻率位置）的模是A3的N/2倍，N為采樣點(diǎn)， 第四個(gè)峰值（頻率位置）的模是A4的N/2倍，N為采樣點(diǎn)， 依次下去...... 考慮到數(shù)量級(jí)較大，一般進(jìn)行歸一化處理，既然第一個(gè)峰值是A1的N倍，那么將每一個(gè)振幅值都除以N即可 FFT具有對(duì)稱性，一般只需要用N的一半，前半部分即可。 4、將振幅譜進(jìn)行歸一化和取半處理 先進(jìn)行歸一化 normalization_y=abs_y/N #歸一化處理（雙邊頻譜）plt.figure()plt.plot(x,normalization_y,'g')plt.title('雙邊頻譜(歸一化)',fontsize=9,color='green')plt.show() 結(jié)果為： 現(xiàn)在我們發(fā)現(xiàn)，振幅譜的數(shù)量級(jí)不大了，變得合理了，接下來進(jìn)行取半處理： half_x = x[range(int(N/2))] #取一半?yún)^(qū)間normalization_half_y = normalization_y[range(int(N/2))] #由于對(duì)稱性，只取一半?yún)^(qū)間（單邊頻譜）plt.figure()plt.plot(half_x,normalization_half_y,'b')plt.title('單邊頻譜(歸一化)',fontsize=9,color='blue')plt.show()這就是我們最終的結(jié)果，需要的“振幅譜”。 三、完整代碼 import numpy as npfrom scipy.fftpack import fft,ifftimport matplotlib.pyplot as pltfrom matplotlib.pylab import mpl mpl.rcParams['font.sans-serif'] = ['SimHei'] #顯示中文mpl.rcParams['axes.unicode_minus']=False #顯示負(fù)號(hào) #采樣點(diǎn)選擇1400個(gè)，因?yàn)樵O(shè)置的信號(hào)頻率分量最高為600赫茲，根據(jù)采樣定理知采樣頻率要大于信號(hào)頻率2倍，所以這里設(shè)置采樣頻率為1400赫茲（即一秒內(nèi)有1400個(gè)采樣點(diǎn)，一樣意思的）x=np.linspace(0,1,1400) #設(shè)置需要采樣的信號(hào)，頻率分量有200，400和600y=7*np.sin(2*np.pi*200*x) + 5*np.sin(2*np.pi*400*x)+3*np.sin(2*np.pi*600*x) fft_y=fft(y) #快速傅里葉變換 N=1400x = np.arange(N) # 頻率個(gè)數(shù)half_x = x[range(int(N/2))] #取一半?yún)^(qū)間 abs_y=np.abs(fft_y) # 取復(fù)數(shù)的絕對(duì)值，即復(fù)數(shù)的模(雙邊頻譜)angle_y=np.angle(fft_y) #取復(fù)數(shù)的角度normalization_y=abs_y/N #歸一化處理（雙邊頻譜） normalization_half_y = normalization_y[range(int(N/2))] #由于對(duì)稱性，只取一半?yún)^(qū)間（單邊頻譜） plt.subplot(231)plt.plot(x,y) plt.title('原始波形') plt.subplot(232)plt.plot(x,fft_y,'black')plt.title('雙邊振幅譜(未求振幅絕對(duì)值)',fontsize=9,color='black') plt.subplot(233)plt.plot(x,abs_y,'r')plt.title('雙邊振幅譜(未歸一化)',fontsize=9,color='red') plt.subplot(234)plt.plot(x,angle_y,'violet')plt.title('雙邊相位譜(未歸一化)',fontsize=9,color='violet') plt.subplot(235)plt.plot(x,normalization_y,'g')plt.title('雙邊振幅譜(歸一化)',fontsize=9,color='green') plt.subplot(236)plt.plot(half_x,normalization_half_y,'blue')plt.title('單邊振幅譜(歸一化)',fontsize=9,color='blue') plt.show() |
|
|
