@zigen 's note
FFT
最終更新:
mynote
FFT(高速フーリエ変換)の使い方
今回使用するFFTライブラリ
y(t) -> Y(ω) : 時間の関数yを周波数空間へ
y(x) -> Y(k):位置の関数yを波数空間へ
Tω=2π, Xk=2πより周波領域、波数領域への変換を簡単に行う。
京都大学の大浦 拓哉先生のFFT
http://www.kurims.kyoto-u.ac.jp/~ooura/fftman/
今回は上記のFFTライブラリの使い方をC言語、Fortranで説明します。
Fortran
データの作成
!======================================================================= ! FFT Test code ! ----- fast fourier transfer test code ----- ! coder : Daisuke Saito ! date : April 11, 2010 ! ! (1)set the initial data about N, pi ! (2)set the initial data for f(i) ! (3)f(i) transfer datas in a(i) ! (4)FFT for a(i) ! (5)out put !======================================================================= program main integer nmax, nmaxsqrt, i parameter (nmax = 32768, nmaxsqrt = 128) integer N, ip(0 : nmaxsqrt + 1) real*8 :: a(0 : nmax), f(0 : nmax), w(0 : nmax*5/4 - 1), pi, dx open(100, file="sin.txt") open(101, file="spectram.txt") write(100,*) "#i, argument, amplitude" write(101,*) "#i, powerspectram" N = nmaxsqrt pi = dacos(-1.d0) dx = 1.d0/N ip(0) = 0 a(:) = 0.d0 do i = 0, N f(i) = dsin(2.d0*pi*i*dx) + dsin(30.d0*2.d0*pi*i*dx) + 3*dsin(7.d0*2.d0*pi*i*dx) !<-この行で初期振動をセット write(100,*) i, 2.d0*pi*i*dx, f(i) enddo do i = 0, N a(2*i) = f(i) !<-Point1 enddo !---Kyoto-u FFT subrotine-- call cdft(2*N, 1, a, ip, w) !<- Point2 do i = 0, N/2 - 1 write(101,*) i, 2.d0*dsqrt(a(2*i)*a(2*i)+a(2*i+1)*a(2*i+1))/N !<-Point3 enddo close(100) close(101) end program main
Point1
今回複素フーリエ変換するのため実数と複素数の入れ物を用意しなければならない、そこでaという入れ物を用意しアーギュメントiが偶数にフーリエ変換したい関数の実数部、iが奇数に虚数部をセットする。
注:今回は虚数部の値がないので、実数部だけ。
Point2
拓哉先生のFFTの説明どおり、cdft(2*N, 1, a, ip, w)を設定
Input: 2*N = 複素FFTをかける配列の大きさ、今回は実部がN個、虚部がN個なので2N個
Input: 1は普通のFFT、-1は逆変換
Input&Output: 変換するデータ、出力データは同じ配列に上書きされる
Input: ?
Output: ワーキング領域
Point3
今回パワースペクトルP(k)を出力したいのでP(k) =|F(k)|^2/N= Fre(k)*Fim(k)/Nを出力
Inputデータ
今回入力するデータは f = sin(x) + sin(30x) + 3sin(7x)
Outputデータ
パワースペクトルのピークがk=1,7,30のところに出ている。
今回使用したデータとか全部
https://docs.google.com/leaf?id=0B08q5BonMYJVMGEzNDkzZGMtOThkMC00OWRkLTk5NzctODVmOTQ2NWUyYmM1&hl=ja