1.題目?

????? 給出時域函數(shù),對其采樣獲得采樣數(shù)據(jù)，并畫出對應(yīng)的x-t曲線；然后對上述采樣數(shù)據(jù)進行傅里葉變換，畫出頻譜圖。運行該程序，需要從鍵盤輸入N中的r值，該程序可以接受的值為5-10的整數(shù)。

---

2.算法及分析

?? 一般地來說，若具有N個采樣點，第L層節(jié)點的值為：

??? 在處理該問題的時候，我們要確定的指數(shù)值。對于節(jié)點,其p值的可以按照下列方式確定。首先將n用r位二進制形式表示出來，再將改二進制右移r-l位，左邊空位補零；然后在將碼序倒置，最后再將改二進制寫出十進制形式，就確定p值。

??? 而上式中涉及的兩個點為對偶節(jié)點。在程序上實現(xiàn)的時候注意對偶節(jié)點。

---

3.程序

import numpy as np import math as ma import cmath import matplotlib.pyplot as plt def f(t,N):x=ma.cos(2*ma.pi/N*t)+0.5*ma.cos(2*2*ma.pi/N*t)+0.8*ma.cos(5*2*ma.pi/N*t)return x r=int(input('輸入（5-10）整數(shù)'));N=2**r;n=N/2;c=0;nu=0; t=np.arange(0,N,1);y=np.zeros((N,r));p=np.zeros((N,r)) x=np.zeros((N,1));x1=[] for i in range(N):m=list(bin(i))for j in range(len(m)-2):y[i,r-1-j]=float(m[len(m)-j-1]) for l in range(1,r+1):z=np.zeros((N,r))for no in range(r):if r-no-1-(r-l)>=0:z[:,r-no-1]=y[:,r-no-1-(r-l)]for nk in range(N):c=0for mk in range(r):c=c+z[nk,mk]*2**mkp[nk,l-1]=c for nk in range(N):x[nk,0]=f(nk,N) x=x+0j for j in range(N):x1.append(j) y1=np.array([x1]); for io in range(1,r+1):y1=y1.reshape((2**io,-1))b=int(y1.shape[0]/2);l=y1.shape[1];lp=0mn=np.zeros((2,int(N/2)))for nk in range(l):for nl in range(b):bg=2*(nl+1)-1bv=2*(nl+1)-2mn[0,lp]=y1[bv,nk]mn[1,lp]=y1[bg,nk]lp=lp+1a1=np.lexsort(mn[::-1,:])mn=mn[:,a1]a2=np.lexsort(mn[0:2,:])mn=mn[:,a2];nk1=np.zeros((N,1));for lop in range(int(N)):nk1[lop,0]=x[lop,0]for nk in range(int(N/2)):pg=[]for gg in mn[:,nk]:pg.append(gg)x[int(pg[0]),0]=nk1[int(pg[0]),0]+cmath.exp(-2*ma.pi*1j/int(N)*p[int(pg[0]),io-1])*nk1[int(pg[1]),0]x[int(pg[1]),0]=nk1[int(pg[0]),0]+cmath.exp(-2*ma.pi*1j/int(N)*p[int(pg[1]),io-1])*nk1[int(pg[1]),0] x=abs(x) bf=x.reshape(1,-1) x=np.arange(0,N,1) x=np.array([x]) y=[];m=0 zk=np.zeros((1,int(N))) for ba in p[:,-1]:zk[0,int(ba)]=bf[0,m]m=m+1 for bb in zk:for bh in bb:y.append(bh)for nm in range(100):y.append(0) vf=[] for nn in x:for vvc in nn:vf.append(vvc)for llk in range(100):vf.append(vvc+llk*0.01) plt.plot(vf,y) plt.rcParams['font.sans-serif']=['SimHei'] plt.rcParams['axes.unicode_minus'] = False plt.xlabel('K值') plt.ylabel('頻譜幅度') plt.title('FFT頻域') plt.show() mk=[];po=[] for bg in range(int(N)+1):mk.append(f(bg,N))po.append(bg) mk1=[];pno=[] for bg2 in range(0,int(N)*25):mk1.append(f(bg2,N))pno.append(bg2) plt.plot(pno,mk1) plt.scatter(po,mk) plt.ylabel('x') plt.xlabel('t') plt.title('FFT時域') plt.show() plt.scatter(po,mk) plt.plot(po,mk) plt.ylabel('x') plt.xlabel('t') plt.title('FFT時域采樣點') plt.show()

---

4.結(jié)果

（1）r=5

---

（2）r=7

---

5.實驗總結(jié)

??? （1）從程序的結(jié)果來看，我們可以用有限的頻譜空間來表征比較復(fù)雜的時域空間，可以將問題大大的簡化；

??? （2）本次程序使用的方法為FFT方法，它的優(yōu)勢很明顯，就是快速。從多次的乘法和加法變?yōu)閹状蔚某朔ê秃唵蔚募訙p法，可以大大的減少計算量；

??? （3）本次程序的編程難點：1.確定p的值2.確定對偶節(jié)點；本次程序較為優(yōu)勢的時候，可以輸出所有P點以及每層的對偶節(jié)點。相對來說更容易展示程序的正確性。

????????????????????????????????????????????

總結(jié)

以上是生活随笔為你收集整理的Python——快速傅里叶变换的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。