INFR

フーリエ変換による積分―INFR
プログラムINFR（Integration by Fouri er Trans form ）は，時刻歴 x ( t ) が，等時間間隔ごとの N個の 離散値データ x m （ m = 0 , 1 , 2 , L , N − 1）として与えられたとき，フーリエ変換およびフーt =m∆tリエ逆 変換を 利用す るこ とによ って， 積分 の離散 値 ∫ x ( t ) dt  （ m = 0 , 1 , 2 , L , N − 1）を0m求めるサブルーチン副プログラムである．

INFR（ フーリエ変換による積分）【目 的】
高速フーリエ変換を利用して，等時間間隔の時刻歴データを積分する．
【使用法】（1）接続方法
CA L LINFR（N， X， ND，DT）

引 数型

NI

プログラムを呼ぶときの内容

プログラムから戻ったときの内容

データの個数 Ｎ≦8192

不 変積分されるデータ

積分したデータRX

1次元配列（ND）ND

I

配列 X の整合寸法

不 変DTR

データの時間間隔

不 変（2）必要なサブルーチンおよび関数副プログラム
F AS T【計算法】
時刻歴 x ( t ) が，等時 間間隔 ∆t ご との N 個 の離散値デ ータ x m （ m = 0 , 1 , 2 , 3 , N − 1）とし て与えられたとき，これをフーリエ変換して，周波数領域における複素フーリエ係数 C k を求め ，次の演算 ( N −1) C 0  π  N / 2−1 J (C k ) 
 −S0 = 22
 Nk k =1
 π   π  i CkS k =  − 1+ i cot   k  C 0   −, S N −k = S k*
 N  Nkπ 
S N / 2 = − C N  0∑を行ったのち-1-
k =1 , 2 , L , N / 2 −1 INFR


∫

N ∆t

x dt  =
m
02πt

N −1

∑S ei ( 2π k m / N )km = 0 , 1 , 2 , L , N −1

k =0によってフーリエ逆変換し，積分の離 散値

tt =m∆t∫ x dt  = ∫ x ( t ) dt
0m0m = 0 , 1 , 2 , L , N −1を求める．
フーリエ変換およびフーリエ 逆変換には，高速フーリエ変 換のプログラム F AS T を使って いる．F AS Tを使うのに都合がよいよう，後続のゼロを付け加えて，総データ数をいったん２の累乗数とし ，また データを 複素数 化する 操作も， プログ ラムの中 の INITIALIZATION のブロ ックで行わ れている ．フ−リ エ変換が 終わった 時点で， 複素フー リエ係数 C k は，す べて N 倍 されていることに，留意が必要である．
【プログラム】
CCC
CCC* * * * * * * * * * * * * * * * * * * * * * * * * * *SUBROUTINE FOR INTEGRATION BY FOURIER TRANSFORM* * * * * * * * * * * * * * * * * * * * * * * * * * *

INFR
INFR
INFR
INFR
CODED BY Y.OHSAKI
INFR
INFR
SUBROUTINE INFR(N,X,ND,DT)
INFR
C
INFR
COMPLEXC(8192)
INFR
DIMENSION X(ND)
INFR
PARAMETER (PI2=6.283185)
INFR
C
INFR
C
INITIALIZATION
INFR
C
INFR
DO 110 M=1,N
INFR
C(M)=CMPLX(X(M),0.)INFR
110 CONTINUE
INFR
NT=2
INFR
120 IF(NT.GE.N) GO TO 130
INFRNT=NT*2
INFR
GO TO 120INFR
130 IF(NT.EQ.N) GO TO 150
INFR
DO 140 M=N+1,NT
INFR
C(M)=(0.,0.)INFR
140 CONTINUE
INFR
150 NFOLD=NT/2+1
INFR
PN=PI2/2./REAL(NT)
INFRC
INFRCFOURIER TRANSFORM
INFR
C
INFR
CALL FAST(NT,C,8192,-1)
INFR
C
INFR
C
INTEGRATION
INFR
C
INFR
C1=REAL(C(1))
INFR
S1=REAL(NT-1)/2.*C1*PN
INFR
DO 160 K=2,NFOLD-1
INFR
S1=S1-AIMAG(C(K))/REAL(K-1)
INFR
C(K)=CMPLX(-1.,1./TAN(REAL(K-1)*PN))*C1*PN-(0.,1.)*C(K)/REAL(K-1) INFR
C(NT-K+2)=CONJG(C(K))
INFR
160 CONTINUE
INFR-2-1234567891011121314151617181920212223242526272829303132333435363738394041INFR

C(1)=CMPLX(S1*2.,0.)C(NFOLD)=CMPLX(-C1*PN,0.)
CCC

INFR
INFR
INFR
INFR
INFR
INFR
INFR
INFR
INFR
INFR
INFR

FOURIER INVERSE TRANSFORMCALL FAST(NT,C,8192,+1)
DO 170 M=1,N
X(M)=REAL(C(M))/PI2*DT
170 CONTINUE
RETURN
END

4243
4445
46
47
48
49505152【使用例 】 加速 度時刻歴 をファイ ル EQ. 01 から読み 取り，サ ブルーチ ン INF R で 積分し て速度時刻歴を求め ，次いでこれを別 掲のサブルーチン DIF R で微分して， もとの加速度時刻 歴に戻せ．
［解］プログラム：DIMENSION DATA(800),VEL(800),ACC(800)
CREAD(5,501) DT,NN,(DATA(M),M=1,NN)
CALL INFR(NN,DATA,800,DT)
DO 110 M=1,NN
VEL(M)=DATA(M)110 CONTINUECALL DIFR(NN,DATA,800,DT)
DO 120 M=1,NN
ACC(M)=DATA(M)120 CONTINUE
STOPC
501 FORMAT(T51,F10.0,I10/(8F10.0))
END

123456789101112131415

アウトプット：積分したデ ータが配列 VEL に，これを微分 したデータが配列 ACC に格納さ れている． これらを，も との加速度時 刻歴，積分し た速度時刻歴 ，微分した 加
速度時刻歴の順序に図示すれば，次図のとおりである．

加速度（gal）

400

最大加速度 326.0 gal

0- 400-3-INFR

速度（kine）

50

最大速度 44.53 kine

0- 50

加速度（gal）

400

最大加速度 325.8 gal

0- 400
0

1

2

3

4
時間（sec）-4-5

678