Re: 数値フーリェ逆変換投稿者：山中和義 投稿日：2014年 2月14日(金)14時38分34秒 |
|
> No.3318[元記事へ] 島村1243さんへのお返事です。 > 解析的にフーリェ変換した複素関数F(ω）を離散値の時間関数f(t)に戻すプログラム(IFFT)をお教えください。 > 例：F(ω)＝1.021*(1/(0.01875+j*ω)-1/(5.01+j*ω)) > 表示時間の全幅：t=0～50[sec] > なお、上記F(ω)の元関数f(t)はf(t)=1.021*(exp(-0.01875*t)-exp(-5.01*t))です。 OPTION ARITHMETIC COMPLEX !複素数モード SET bitmap SIZE 960,160 !表示領域 SET WINDOW -2,52,-1,2 DRAW grid(2,0.5) LET T=50 !時間区間 [0,T] DEF G(jw)=1.021*(1/(0.01875+jw)-1/(5.01+jw)) LET N=1024*8 !データ総数 ※大きいほど精度がよい DIM d(0 TO N-1) !入出力用配列 LET Gamma=5 !※|γ|=3～7 !γを大きくすれば精度は上がりそうだが、 !あとでExp(γt)をかけるので、tが大きいところで発散する。 LET r=Gamma/T !γを決める FOR k=0 TO N/2 !変換データの作成 LET Fs=G( COMPLEX(r, 2*PI*k/T) ) !sk=γ+i*2π*k/T、k=0～n-1 LET Hw=(COS(2*PI*k/N)+1)/2 !データは離散なため、ハニング関数をかけて平滑化する ※精度の向上 LET d(k)=N/T*Fs*Hw !n/T*Y(k) * Hw 0～N/2 IF NOT(k=0 OR k=N/2) THEN LET d(N-k)=COMPLEX(Re(d(k)),-Im(d(k))) !N-1～N/2+1は、共役複素数 NEXT k CALL IFFT(d) !高速逆フーリエ変換 LET dt=T/N !時間刻み幅 ⊿t FOR k=0 TO N-1 STEP 8 !結果の表示 [0,T] LET tt=k*dt !経過時間 PLOT LINES: tt, Re(d(k))*EXP(r*tt); !Exp(γt)をかける ※Exp(γ*k*T/n)、k=0～n-1 NEXT k PLOT LINES !検算 DEF f(t)=1.021*(EXP(-0.01875*t)-EXP(-5.01*t)) SET LINE COLOR 4 FOR k=0 TO T STEP 1/2^8 PLOT LINES: k,f(k); NEXT k PLOT LINES END EXTERNAL SUB IFFT(x()) !高速逆フーリエ変換 x() : 入力/出力データ OPTION ARITHMETIC COMPLEX !複素数モード DECLARE EXTERNAL SUB FFTMAIN LET nx=SIZE(x) LET theta=2*PI/nx CALL FFTMAIN(x, theta) MAT x=(1/nx)*x END SUB EXTERNAL SUB FFTMAIN(x(), theta) OPTION ARITHMETIC COMPLEX !複素数モード LET nx=SIZE(x) IF MOD(nx, 2)<>0 THEN !DFTの計算 DIM w(0 TO nx-1), xtmp(0 TO nx-1) MAT xtmp=x FOR k=0 TO nx-1 LET tmp=theta*k FOR n=0 TO nx-1 LET w(n)=EXP( COMPLEX(0, -tmp*n) ) NEXT n LET x(k)=DOT(w, xtmp) NEXT k ELSE !２分して再帰呼出し LET hnx=nx/2 DIM x0(0 TO hnx-1), x1(0 TO hnx-1) FOR k=0 TO hnx-1 LET x0(k)=x(k)+x(k+hnx) LET wk=EXP( COMPLEX(0, theta*k) ) LET x1(k)=wk*(x(k)-x(k+hnx)) NEXT k CALL FFTMAIN(x0, 2*theta) CALL FFTMAIN(x1, 2*theta) FOR k=0 TO hnx-1 LET x(2*k)=x0(k) LET x(2*k+1)=x1(k) NEXT k END IF END SUB 
|