1.题目
给出时域函数,对其采样获得采样数据,并画出对应的x-t曲线;然后对上述采样数据进行傅里叶变换,画出频谱图。运行该程序,需要从键盘输入N中的r值,该程序可以接受的值为5-10的整数。
2.算法及分析
一般地来说,若具有N个采样点,第L层节点的值为:
在处理该问题的时候,我们要确定的指数值。对于节点,其p值的可以按照下列方式确定。首先将n用r位二进制形式表示出来,再将改二进制右移r-l位,左边空位补零;然后在将码序倒置,最后再将改二进制写出十进制形式,就确定p值。
而上式中涉及的两个点为对偶节点。在程序上实现的时候注意对偶节点。
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)整数));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=0 for mk in range(r): c=c+z[nk,mk]*2**mk p[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=0 mn=np.zeros((2,int(N/2))) for nk in range(l): for nl in range(b): bg=2*(nl+1)-1 bv=2*(nl+1)-2 mn[0,lp]=y1[bv,nk] mn[1,lp]=y1[bg,nk] lp=lp+1 a1=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.结果
(1)r=5
(2)r=7
5.实验总结
(1)从程序的结果来看,我们可以用有限的频谱空间来表征比较复杂的时域空间,可以将问题大大的简化;
(2)本次程序使用的方法为FFT方法,它的优势很明显,就是快速。从多次的乘法和加法变为几次的乘法和简单的加减法,可以大大的减少计算量;
(3)本次程序的编程难点:1.确定p的值2.确定对偶节点;本次程序较为优势的时候,可以输出所有P点以及每层的对偶节点。相对来说更容易展示程序的正确性。