本文實例為大家分享了python傅里葉變換FFT繪制頻譜圖的具體代碼,供大家參考,具體內(nèi)容如下
創(chuàng)新互聯(lián)是一家專業(yè)從事成都做網(wǎng)站、成都網(wǎng)站設(shè)計、網(wǎng)頁設(shè)計的品牌網(wǎng)絡(luò)公司。如今是成都地區(qū)具影響力的網(wǎng)站設(shè)計公司,作為專業(yè)的成都網(wǎng)站建設(shè)公司,創(chuàng)新互聯(lián)依托強大的技術(shù)實力、以及多年的網(wǎng)站運營經(jīng)驗,為您提供專業(yè)的成都網(wǎng)站建設(shè)、營銷型網(wǎng)站建設(shè)及網(wǎng)站設(shè)計開發(fā)服務(wù)!頻譜圖的橫軸表示的是 頻率, 縱軸表示的是振幅
#coding=gbk import numpy as np import pandas as pd import matplotlib.pyplot as plt #依據(jù)快速傅里葉算法得到信號的頻域 def test_fft(): sampling_rate = 8000 #采樣率 fft_size = 8000 #FFT長度 t = np.arange(0, 1.0, 1.0/sampling_rate) x = np.sin(2*np.pi*156.25*t) + 2*np.sin(2*np.pi*234.375*t)+ 3*np.sin(2*np.pi*200*t) xs = x[:fft_size] xf = np.fft.rfft(xs) / fft_size #返回fft_size/2+1 個頻率 freqs = np.linspace(0, sampling_rate/2, fft_size/2+1) #表示頻率 xfp = np.abs(xf) * 2 #代表信號的幅值,即振幅 plt.figure(num='original', figsize=(15, 6)) plt.plot(x[:100]) plt.figure(figsize=(8,4)) plt.subplot(211) plt.plot(t[:fft_size], xs) plt.xlabel(u"時間(秒)", fontproperties='FangSong') plt.title(u"156.25Hz和234.375Hz的波形和頻譜", fontproperties='FangSong') plt.subplot(212) plt.plot(freqs, xfp) plt.xlabel(u"頻率(Hz)", fontproperties='FangSong') plt.ylabel(u'幅值', fontproperties='FangSong') plt.subplots_adjust(hspace=0.4) plt.show() test_fft() # np.clip(a, a_min, a_max, out) 輸出與a 的shape一樣,大于等于a_min,小于等于a_max的數(shù),即在 [a_min, a_max]之間的數(shù) a = np.arange(10) print(a) print(a.shape) # [0 1 2 3 4 5 6 7 8 9] b = np.empty((10,)) np.clip(a, 3, 8, out=b) print(b) # [3. 3. 3. 3. 4. 5. 6. 7. 8. 8.] c = np.clip(a, 4, 10) print(c) # [4 4 4 4 4 5 6 7 8 9] #a_min, a_max也可以輸入與a 相同shape的數(shù)組 d = np.arange(4) d1 = np.clip(d, [-1, 1, -3, 2], 2) print(d) print(d1) # [0 1 2 3] #原數(shù)組 # [0 1 2 2] print(np.log10(1000)) def test_fft(): # FFT變換是針對一組數(shù)值進行運算的,這組數(shù)的長度N必須是2的整數(shù)次冪,例如64, 128, 256等等; 數(shù)值可以是實數(shù)也可以是復(fù)數(shù), # 通常我們的時域信號都是實數(shù),因此下面都以實數(shù)為例。我們可以把這一組實數(shù)想像成對某個連續(xù)信號按照一定取樣周期進行取樣而得來, # 如果對這組N個實數(shù)值進行FFT變換,將得到一個有N個復(fù)數(shù)的數(shù)組,我們稱此復(fù)數(shù)數(shù)組為頻域信號,此復(fù)數(shù)數(shù)組符合如下規(guī)律: # # 下標(biāo)為0和N/2的兩個復(fù)數(shù)的虛數(shù)部分為0, # 下標(biāo)為i和N-i的兩個復(fù)數(shù)共軛,也就是其虛數(shù)部分?jǐn)?shù)值相同、符號相反。 np.random.seed(66) X = np.random.rand(8) print(X) # [0.15428758 0.13369956 0.36268547 0.67910888 0.19445006 0.25121038 # 0.75841639 0.55761859] xf = np.fft.fft(X) print(xf) # [ 3.0914769 +0.j -0.20916178+0.39291702j -0.77236422+0.85181752j # 0.12883683-0.39854483j -0.15179792+0.j 0.12883683+0.39854483j # -0.77236422-0.85181752j -0.20916178-0.39291702j] #通過快速傅里葉變換的逆變換 ifft 還原成原來的值 X1 = np.fft.ifft(xf) print(X1) # [0.15428758+0.00000000e+00j 0.13369956-2.00387919e-16j # 0.36268547+1.66533454e-16j 0.67910888+1.51815661e-16j # 0.19445006+0.00000000e+00j 0.25121038-1.51815661e-16j # 0.75841639-1.66533454e-16j 0.55761859+2.00387919e-16j] # 下面讓我們來看看FFT變換之后的那些復(fù)數(shù)都代表什么意思。 # # 首先下標(biāo)為0的實數(shù)表示了時域信號中的直流成分的多少 # 下標(biāo)為i的復(fù)數(shù)a+b*j表示時域信號中周期為N/i個取樣值的正弦波和余弦波的成分的多少, 其中a表示cos波形的成分,b表示sin波形的成分 X = np.ones(8) x2 = np.fft.fft(X) / len(X) # 為了計算各個成分的能量多少,需要將FFT的結(jié)果除以FFT的長度 print(x2) # [1.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j 0.+0.j] X = np.arange(0, 2*np.pi, 2*np.pi/8) y = np.sin(X) x3 = np.fft.fft(y) /len(y) print(x3) # [ 1.43029718e-18+0.00000000e+00j -4.44089210e-16-5.00000000e-01j # 只有下標(biāo)為 1 的復(fù)數(shù)的虛部為-0.5, # 1.53080850e-17-1.38777878e-17j 3.87727691e-17-1.11022302e-16j # 2.91858728e-17+0.00000000e+00j 0.00000000e+00-1.11022302e-16j # 1.53080850e-17+1.38777878e-17j 3.44084101e-16+5.00000000e-01j] output1 = np.fft.fft(np.cos(X) / len(X)) print(output1) # [-4.30636606e-17+0.00000000e+00j 5.00000000e-01-2.66538563e-16j #只有下標(biāo)為1 的實部為 0.5 # 1.53080850e-17+0.00000000e+00j 5.55111512e-17+1.97149624e-16j # 1.24474906e-17+0.00000000e+00j -1.11022302e-16+2.05306223e-16j # 1.53080850e-17+0.00000000e+00j 5.00000000e-01-1.35917284e-16j] #綜合的例子 X = np.arange(0, 2*np.pi, 2*np.pi/128) y = 0.3*np.cos(X) + 0.5*np.cos(2*X+np.pi/4) + 0.8*np.cos(3*X-np.pi/3) yf = np.fft.fft(y) / len(y) print(2*np.abs(yf[1]), np.rad2deg(np.angle(yf[1]))) # 0.30000000000000016 3.3130777931911615e-15 #計算出幅值和相位角 print(2*np.abs(yf[2]), np.rad2deg(np.angle(yf[2]))) # 0.5000000000000002 44.999999999999986 print(2*np.abs(yf[3]), np.rad2deg(np.angle(yf[3]))) # 0.7999999999999998 -60.00000000000007 # 周期為128/1.0點的余弦波的相位為0, 振幅為0.3 # 周期為64/2.0點的余弦波的相位為45度, 振幅為0.5 # 周期為128/3.0點的余弦波的相位為-60度,振幅為0.8 # test_fft() #使用多個正玄波合成三角波 import pylab as pl # 取FFT計算的結(jié)果freqs中的前n項進行合成,返回合成結(jié)果,計算loops個周期的波形 def fft_combine(freqs, n, loops=1): length = len(freqs) * loops data = np.zeros(length) index = loops * np.arange(0, length, 1.0) / length * (2 * np.pi) for k, p in enumerate(freqs[:n]): if k != 0: p *= 2 # 除去直流成分之外,其余的系數(shù)都*2 data += np.real(p) * np.cos(k*index) # 余弦成分的系數(shù)為實數(shù)部 data -= np.imag(p) * np.sin(k*index) # 正弦成分的系數(shù)為負(fù)的虛數(shù)部 return index, data # 產(chǎn)生size點取樣的三角波,其周期為1 def triangle_wave(size): x = np.arange(0, 1, 1.0/size) y = np.where(x<0.5, x, 0) y = np.where(x>=0.5, 1-x, y) return x, y def test_show(): fft_size = 256 # 計算三角波和其FFT x, y = triangle_wave(fft_size) fy = np.fft.fft(y) / fft_size # 繪制三角波的FFT的前20項的振幅,由于不含下標(biāo)為偶數(shù)的值均為0, 因此取 # log之后無窮小,無法繪圖,用np.clip函數(shù)設(shè)置數(shù)組值的上下限,保證繪圖正確 pl.figure() pl.plot(np.clip(20*np.log10(np.abs(fy[:20])), -120, 120), "o") pl.xlabel("frequency bin") pl.ylabel("power(dB)") pl.title("FFT result of triangle wave") # 繪制原始的三角波和用正弦波逐級合成的結(jié)果,使用取樣點為x軸坐標(biāo) pl.figure() pl.plot(y, label="original triangle", linewidth=2) for i in [0,1,3,5,7,9]: index, data = fft_combine(fy, i+1, 2) # 計算兩個周期的合成波形 pl.plot(data, label = "N=%s" % i) pl.legend() pl.title("partial Fourier series of triangle wave") pl.show() # test_show()
另外有需要云服務(wù)器可以了解下創(chuàng)新互聯(lián)scvps.cn,海內(nèi)外云服務(wù)器15元起步,三天無理由+7*72小時售后在線,公司持有idc許可證,提供“云服務(wù)器、裸金屬服務(wù)器、高防服務(wù)器、香港服務(wù)器、美國服務(wù)器、虛擬主機、免備案服務(wù)器”等云主機租用服務(wù)以及企業(yè)上云的綜合解決方案,具有“安全穩(wěn)定、簡單易用、服務(wù)可用性高、性價比高”等特點與優(yōu)勢,專為企業(yè)上云打造定制,能夠滿足用戶豐富、多元化的應(yīng)用場景需求。
當(dāng)前文章:python傅里葉變換FFT繪制頻譜圖-創(chuàng)新互聯(lián)
URL網(wǎng)址:http://www.rwnh.cn/article44/dghihe.html
成都網(wǎng)站建設(shè)公司_創(chuàng)新互聯(lián),為您提供網(wǎng)頁設(shè)計公司、自適應(yīng)網(wǎng)站、關(guān)鍵詞優(yōu)化、網(wǎng)站排名、企業(yè)建站、手機網(wǎng)站建設(shè)
聲明:本網(wǎng)站發(fā)布的內(nèi)容(圖片、視頻和文字)以用戶投稿、用戶轉(zhuǎn)載內(nèi)容為主,如果涉及侵權(quán)請盡快告知,我們將會在第一時間刪除。文章觀點不代表本網(wǎng)站立場,如需處理請聯(lián)系客服。電話:028-86922220;郵箱:631063699@qq.com。內(nèi)容未經(jīng)允許不得轉(zhuǎn)載,或轉(zhuǎn)載時需注明來源: 創(chuàng)新互聯(lián)
猜你還喜歡下面的內(nèi)容