主要用于地震動(dòng)處理種的代碼：

傅里葉譜 def myfft(inp,delta): #輸入的x為豎向的arrayfrom scipy.fftpack import fft,ifftinp=np.array(inp).reshape(-1,1)inpf=np.fft.fft(inp,axis=0) #默認(rèn)的是對(duì)行向量fourier變換，所以此處要定義下軸inpf=inpf*2/inp.shape[0]inpf2=inpf[0:int((inp.shape[0]+1)/2)]#取一半n_points=inpf2.shape[0]#地震波點(diǎn)數(shù)frequency=1/deltafseries=np.linspace(0,frequency/2,n_points)# plt.plot(fseries,inpf2)return inpf2,fseries

?

功率譜 def mypowverspectra(inp,delta):from scipy.fftpack import fft,ifftinp=np.array(inp).reshape(-1,1)inpf=np.fft.fft(inp,axis=0) #默認(rèn)的是對(duì)行向量fourier變換，所以此處要定義下軸inpf=inpf*2/inp.shape[0]inpf2=inpf[0:int((inp.shape[0]+1)/2)]#取一半inpf2=np.abs(inpf2)**2n_points=inpf2.shape[0]#地震波點(diǎn)數(shù)frequency=1/deltafseries=np.linspace(0,frequency/2,n_points)# plt.plot(fseries,inpf2)return inpf2,fseries

?

反應(yīng)譜 def myresponse(inpacc,delta,damp,Tmax):#Tmax是反應(yīng)譜周期最大值inpacc=np.array(inpacc).reshape(-1,1)count=inpacc.shape[0]displace=np.zeros([count])velocity=np.zeros([count])absacce=np.zeros([count])ta=np.arange(delta,Tmax,delta) # 頻段范圍-1/deltaHZmdis=np.zeros([len(ta)])#只記錄一種阻尼比的響應(yīng)幅值mvel=np.zeros([len(ta)])macc=np.zeros([len(ta)])frcy=2*math.pi/tadamfrcy=frcy*np.sqrt(1-damp**2)e_t=np.exp(-damp*frcy*delta)s=np.sin(damfrcy*delta)c=np.cos(damfrcy*delta)d_f=(2*damp**2-1)/(frcy**2*delta)d_3t=damp/(frcy**3*delta)for i in range(len(ta)):A=np.zeros([2,2])A[0,0]=e_t[i]*(s[i]*damp/np.sqrt(1-damp**2)+c[i])A[0,1]=e_t[i]*s[i]/damfrcy[i]A[1,0]=-frcy[i]*e_t[i]*s[i]/np.sqrt(1-damp**2)A[1,1]=e_t[i]*(-s[i]*damp/np.sqrt(1-damp**2)+c[i])B=np.zeros([2,2])B[0,0]=e_t[i]*((d_f[i]+damp/frcy[i])*s[i]/damfrcy[i]+(2*d_3t[i]+1/frcy[i]**2)*c[i])-2*d_3t[i]B[0,1]=-e_t[i]*(d_f[i]*s[i]/damfrcy[i]+2*d_3t[i]*c[i])-1/frcy[i]**2+2*d_3t[i]B[1,0]=e_t[i]*((d_f[i]+damp/frcy[i])*(c[i]-damp/np.sqrt(1-damp**2)*s[i])-(2*d_3t[i]+1/frcy[i]**2)*(damfrcy[i]*s[i]+damp*frcy[i]*c[i]))+1/(frcy[i]**2*delta)B[1,1]=e_t[i]*(1/(frcy[i]**2*delta)*c[i]+s[i]*damp/(frcy[i]*damfrcy[i]*delta))-1/(frcy[i]**2*delta)for k in range(0,count-1):displace[k+1]=A[0,0]*displace[k]+A[0,1]*velocity[k]+B[0,0]*inpacc[k]+B[0,1]*inpacc[k+1]velocity[k+1]=A[1,0]*displace[k]+A[1,1]*velocity[k]+B[1,0]*inpacc[k]+B[1,1]*inpacc[k+1]absacce[k+1]=-2*damp*frcy[i]*velocity[k+1]-frcy[i]**2*displace[k+1]mdis[i]=np.max(np.abs(displace))mvel[i]=np.max(np.abs(velocity))if i==0:macc[i]=np.max(np.abs(inpacc))else:macc[i]=np.max(np.abs(absacce)) # plt.plot(ta,macc)return macc,ta

反應(yīng)譜是matlab代碼修改的算法，詳細(xì)見http://blog.sciencenet.cn/blog-419879-446739.html

總結(jié)

以上是生活随笔為你收集整理的python程序精确法求解反应谱,傅里叶谱，功率谱的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

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