離散フーリエ変換と高速フーリエ変換(fortran90)

数値計算で離散フーリエ変換を行う方法です。
言語は
fortran90

intel®のMKL(マス・カーネル・ライブラリー)
を用いて離散フーリエ変換を行いたいと考えます。


高速フーリエ変換は離散フーリエ変換を高速に行う手法のことです。

通常の離散フーリエの計算量のオーダーは、\(O(N^2)\)であり、
基数2の高速フーリエの計算量のオーダーは、\(O(Nlog_2 N)\)です。

基数2は良く言われる通常の配列の個数\(2^m\)個のことを表しています。

離散フーリエ変換を行う際には周期境界条件が課されます。
分点の端で周期境界条件が満たされていない場合、それはその点でステップ関数が含まれることを意味します。
依って予期しない高周波成分が出現するので、注意してください。
一番良い判定方法は関数の対数微分が一致するかどうかを判定するのが良いでしょう。

MKLの離散フーリエ変換ルーチンの概要1


MKLの離散フーリエ変換ルーチンは非常に優秀です。
渡す配列サイズによってプログラムが自動的に判断し、最適な手法で離散フーリエ変換を行います。
MKLのマニュアルによると、

注 : DFT関数は任意の長さをサポートしている。
これらのルーチンは、基数 2 だけでなく、基数 3、5、7、11 に対しても高
性能で広範な機能性を提供する。対応する基数の一覧は、使用しているラ
イブラリー・バージョンの「インテルMKLテクニカル・ユーザー・ノー
ト」を参照のこと。

とあります。
高速フーリエ変換の考え方は\(2^m\)の時だけでなく、\(3^m, 5^m, 7^m, 11^m, \cdots\)の時にも同様に考えることができます。
この時、何の整数乗にするかを区別するために “基数” という言葉を使います。
プログラムを作る場合、通常は一番シンプルな基数”2″が選ばれます。

具体的に配列サイズ\(N\)が
\(N= 2^m\)、ただしmは正の整数(配列をAとするとサイズ\(A(0:2^m-1)\))の時に基数2の高速フーリエ変換が行われます。
それ以外の時は通常の離散フーリエ変換が行われているはずです。

余り自信が無いのではっきりとしたことは言えませんが、分点量(2~20000程度まで)に対するMKLのフーリエ変換の計算速度はほとんどわからないくらいでした。特定の基数の時だけ早いんだろうと予想していましたが、そんなことはなく、それだけ優秀であるという事でしょう。

もう6年前に書かれている記事ですので今はどうなっているかわかりませんが、FFTWと計算速度を比較したページがありましたので載せておきます。
一言でいえば、MKLによる離散フーリエ変換はFFTWよりも早いよ、ということです。
インテル謹製の数値演算ライブラリ「MKL」を使ってプログラムを高速化 3ページ -OSDN, (2009)

MKLの離散フーリエ変換ルーチンの概要2


MKLでの1次元、長さnの順方向離散フーリエ変換は、
\(
\displaystyle z_k=\sum_{j=0}^{n-1}w_je^{-i2\pi jk/n}
\)

に従い行われます(MKLリファレンスマニュアル、11-26)。
数値計算の世界では、余分な係数\(\frac{1}{2\pi}\)をなくすため、指数関数の肩に\(2\pi\)が掛けられたものが良く使用されます。

※ここでは、フーリエ変換前の空間を”位置空間“と呼び、フーリエ変換後の空間を”波数空間“と呼ぶことにします。
位置波数ではなく、時間周波数としても同じことなのでどちらでもいいのですが、ここでは、位置波数の呼び方に統一いたします。

DFTexampleNeven
上図のように位置空間上で区間\(L\)の領域で\(N\)個に分割します。
上図はNが偶数の場合です。
MKLの離散フーリエ変換ルーチンはいかなるNに対してもフーリエ変換が可能ですが、奇数の場合、以下のように波数成分の格納が変わります。
周波数格納順

\(
f(x)=\sin(\alpha x)
\)

を離散フーリエ変換した場合、
波数空間でのフーリエ変換後のピークは\(k=\frac{\alpha}{2\pi}\)に現れます。
この時、f(x)のフーリエ変換後の関数\(g(k)=F_{[f(x)]}\)は複素数として格納され、上記の場合では
\(
\displaystyle g(k)=i\left\{\frac{N}{2}\delta(k-\frac{\alpha}{2\pi})+\frac{N}{2}\delta(k+\frac{\alpha}{2\pi})\right\}
\)

となります。ここで\(i\)は虚数単位を表し、\(\delta(k)\)はデルタ関数を表します。
波数空間での最大の周波数値\(f_c\)はナイキストの定理より、
\(
\begin{align}
f_c = \frac{1}{2\Delta}, \Delta=\frac{L}{N}
\end{align}
\)
と表されます。これは、波数空間では\(\pm f_c\)までサンプリング可能であることを示しています。
フーリエ変換後の値\(g(k)\)の分点上の和\(S\)、すなわち
\(
\displaystyle S=\sum_{j=0}^{N-1}g(k_j)
\)

は離散フーリエ変換の性質により、\(S\)は分点の総数\(N\)に等しくなります。

逆離散フーリエ変換によって元の空間に戻る場合、規格化因子\(N\)で割らなければなりません。

MKLの1次元離散フーリエ変換ルーチンのプログラム例


使い方はFFT Code Examples -Intel®Developer Zone を参考としました。

例として
関数
\(
f(x)=\sin(4x)
\)
をフーリエ変換し、元に戻すプログラムを書きましょう。
プログラムの流れとしては、

関数\(f(x)\)の位置空間での書き出し(ファイル”before.d”に書き出し、横軸:位置\(x\), 縦軸\(f(x)\))
↓ 順方向フーリエ変換(規格化因子で割られない)
関数\(g(k)=F_{[f(x)]}\)の書き出し(ファイル”k.d”に書き出し、横軸:\(k/(2\pi)\), 縦軸\(g(k)\))
↓ 逆方向フーリエ変換(規格化因子Nで割られる)
関数\(f(x)=F^{-1}_{[g(k)]}\)の書き出し(ファイル”after.d”に書き出し、横軸:x, 縦軸f(x))
になっています。

コンパイルを行うにあたり重要な事が2点あります。
1, MKLを使う事
2, program main..のにinclude文を付け加える事(かも?)

このincludeは、MKLにある離散フーリエ変換ルーチンを使いますよーというサインであり、パスはMKLのバージョンや、インストールした状況によって変わります。
うまく、ファイル”mkl_dfti.f90″を見つけてください。パスを上から辿っていくのがいいと思います。

※2016/02/27時点での最新バージョン(ifort –version で確認できます)
ifort (IFORT) 16.0.2 20160204
において、mkl_dfti.fのデフォルトでの場所は、
“/opt/intel/compilers_and_libraries_2016/linux/mkl/include/mkl_dfti.f90”
です。また、このバージョンだけかどうかわかりませんが、include文が無くてもどうやら動くようです。

プログラムのコンパイルと実行は

$ ifort -mkl main.f90
$ ./a.out

で良いでしょう。

include "/opt/intel/mkl/include/mkl_dfti.f90"
program main
  implicit none
  integer::i,n
  double precision,allocatable::x(:),k(:),mk(:)
  complex(kind(0d0)),allocatable::z(:),mz(:)
  double precision::hx,xmin,xmax
 
  double precision,parameter::pi=dacos(-1.d0)
  complex(kind(0d0))::func
  external::func
 
  n=130
  allocate(x(0:n-1),k(0:n-1),z(0:n-1))
  x=0d0; k=0d0; z=dcmplx(0d0,0d0)

  allocate(mk(0:n-1),mz(0:n-1))
  mk=0d0; mz=dcmplx(0d0,0d0)
 
  xmin=-3d0*pi
  xmax=3d0*pi
  hx=(xmax-xmin)/dble(n)

  do i=0,n-1
     x(i)=xmin+dble(i)*hx
     z(i)=func(x(i))
  enddo
 
  open(21,file="before.d")
  do i=0,n-1
     write(21,'(3e20.10e3)')x(i),dble(z(i)),dimag(z(i))
  enddo
  close(21)
 
  call dftf(n,k,hx) ! Get frequency k
 
  call dft(n,z,"forward") ! forward dft

  call dfts(n,k,z,mk,mz) !sort frequency
  open(21,file="k.d")
  do i=0,n-1
     write(21,'(3e20.10e3)')mk(i),dble(mz(i)),dimag(mz(i))
  enddo
  close(21)
 
  call dft(n,z,"backward") ! backward dft
  open(21,file="after.d")
  do i=0,n-1
     write(21,'(3e20.10e3)')x(i),dble(z(i)),dimag(z(i))
  enddo
  close(21)    
 
  stop
end program main

function func(x)
  implicit none
  double precision::x
  complex(kind(0d0))::func
 
  func=dcmplx(dsin(4d0*x),0d0)
 
  return
end function func

上で用いているサブルーチン
dft,dfts,dftfはそれぞれフーリエ変換、フーリエ変換後の空間の値、ソート用のルーチンであることを意味し、その中身は以下のようなものです。例のごとく、MKLのルーチンはそのままでは使いにくいので1クッション挟みます。
サブルーチンのコンセプトは、
ただ順/逆方向離散フーリエ変換をしてくれるサブルーチンです。
それ以外の要素は省いています。各ルーチンの詳細は下に書きました。

subroutine dft(n,z,FB)
  !sikinote
  ! developer : sikino
  ! date : 2015/09/20
  use MKL_DFTI
  implicit none
  integer,intent(in)::n
  complex(kind(0d0)),intent(inout)::z(0:n-1)  
  character(*)::FB
 
  integer::Status
  TYPE(DFTI_DESCRIPTOR),POINTER::hand
 
  !DFT : Discrete Fourier Transform
  !
  !n    --> number of data.
  !z(i) --> value of data at x(i)
  !FB   --> "forward"  : Forward DFT
  !     --> "backward" : Backward DFT
 
  Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,n)
  Status = DftiCommitDescriptor(hand)
  if(trim(FB).eq."forward")then
     Status = DftiComputeForward(hand,z)
  elseif(trim(FB).eq."backward")then
     Status = DftiComputeBackward(hand,z)
  else
     write(6,*)"DFT string different"
     stop
  endif
  Status = DftiFreeDescriptor(hand)

  if(trim(FB).eq."backward")then
      z=z/dble(n)
  end if
       
  return
end subroutine dft

subroutine dfts(N,f,z,mf,mz)
  !sikinote
  ! developer : sikino
  ! date : 2015/09/20
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::f(0:N-1)
  double precision,intent(out)::mf(0:N-1)
  complex(kind(0d0)),intent(in)::z(0:N-1)
  complex(kind(0d0)),intent(out)::mz(0:N-1)

  integer::i,j
 
  if(mod(n,2).eq.0)then
     do i=0,N-1
        if(i.le.N/2)then
           j=i+N/2-1
           mf(j)=f(i)
           mz(j)=z(i)
        else
           j=i-N/2-1
           mf(j)=f(i)
           mz(j)=z(i)
        endif
     enddo
  else
     do i=0,N-1
        if(i.le.(N-1)/2)then
           j=i+(N-1)/2
           mf(j)=f(i)
           mz(j)=z(i)
        else
           j=i-(N-1)/2-1
           mf(j)=f(i)
           mz(j)=z(i)
        endif
     enddo
  endif

  return
end subroutine dfts

subroutine dftf(N,f,h)
  !sikinote
  ! developer : sikino
  ! date : 2015/09/20
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::h
  double precision,intent(out)::f(0:N-1)

  integer::i
  double precision::mf(0:N-1)
 
  if(mod(N,2).eq.0)then
     do i=0,N-1
        mf(i)=(dble(i+1)-dble(N)*0.5d0)/(dble(N)*h)
     enddo

     do i=0,N-1
        if(i.le.N/2-2)then
           f(i+N/2+1)=mf(i)
        else
           f(i-N/2+1)=mf(i)
        endif
     enddo
  else
     do i=0,N-1
        mf(i)=(dble(i)-dble(N-1)*0.5d0)/(dble(N)*h)
     enddo

     do i=0,N-1
        if(i.le.(N-1)/2-1)then
           f(i+(N-1)/2+1)=mf(i)
        else
           f(i-(N-1)/2)=mf(i)
        endif
     enddo
  endif
 
  return
end subroutine dftf

ここ↑に書いてあるサブルーチン

  1. dft(N,z,FB)
  2. dftf(N,f,h)
  3. dfts(N,f,z,mf,mz)

を説明します。

  1. dft(N,z,FB)

    このルーチンによって順/逆方向フーリエ変換が行われます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (入出力) z : complex(kind(0d0)), データ配列z(0:N-1)
    • (入力) FB : character, 順方向(“forward”)、逆方向(“backward”)の指定

    です。データ配列に上書きして値が返ってきます。
    サイズNのデータ配列z,順方向の離散フーリエ変換をしたい場合は、

    call dft(size(z,1),z,"forward")

    とするのがいいでしょう。
    規格化は順方向の時には行わず、逆方向の時にだけ\(1/N\)を掛けています。

  2. dftf(N,f,h)

    このルーチンは波数空間での横軸の値f(0:N-1)を与えます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (出力) f : double precision, 波数空間の横軸の値f(0:N-1)
    • (入力) h : double precision, 位置空間でのサンプルレート(刻み幅)

    です。
    波数fの順番は少しややこしく、
    周波数格納順
    となっています。
    また、ここで出力されるfは波数を\(2\pi\)で割ったものが出力されます。
    具体的には、\(\sin(x)\)をMKLのフーリエ変換を行いグラフを書くとフーリエ変換後の空間(波数空間)で、ピークの位置は\(\pm 0.159115494309\cdots (\sim \frac{1}{2\pi})\)に表れることを意味します。
    また、変換先の空間をグラフで描画したい時に値を小さい順に並べ替えたい必要が出てきたときは並べ替え用のルーチンdftsを使用してください。

  3. dfts(N,f,z,mf,mz)

    このルーチンは変換先の波数空間でグラフを描画したい時に、順番を値を小さい順に並べ替えます。
    引数はそれぞれ

    • (入力) N : integer, データ配列のサイズ
    • (入力) f : double precision, 波数空間の横軸の値f(0:N-1)
    • (入力) z : complex(kind(0d0)), 波数空間でのデータ配列z(0:N-1)
    • (出力) mf : double precision, 並び替えられた波数空間の横軸の値mf(0:N-1)
    • (出力) mz : complex(kind(0d0)), 並び替えられた波数空間でのデータ配列z(0:N-1)

    です。上書きして引数を返す必要性がないだろう、と判断して別の配列に書き出させています。
    もしもメモリが足らないとか、気に食わない人は書き換えてください。


スポンサーリンク


出力ファイル”before.d”, “k.d”, “after.d”を書いてみます。
位置空間での区間\([-3\pi\sim3\pi]\)で、分点数\(N=288\)で行った場合、
位置空間での波形、この場合は\(\sin(4x)\)が描かれます。
位置空間

波数空間では
\(
\displaystyle g(k)=i\left\{\frac{288}{2}\delta(k-\frac{4}{2\pi})+\frac{288}{2}\delta(k+\frac{4}{2\pi})\right\}
\)

が得られるはずで、実際にプロットして確かめてみると
N288_c
となります。
※Nは2以上のどんな整数でもokです。偶数、奇数、素数であっても。

※実部に現れる小さなピークに関して
原因はよくわかりません。\(N=2^m\)ではない時になんかこれが出てきます。
不安な人は\(N=2^m\)に固定して使ってください。
N=1024で行うとこうなります。
N1024

ベンチマーク

同じ計算機上でのベンチマーク結果を載せます。
コンパイラ:ifort (IFORT) 16.0.2 20160204
MKL:/opt/intel/compilers_and_libraries_2016.2.181
計算スレッド数1
ベンチマークdft_mkl_c

ベンチマーク用プログラムはこちら↓

2次元のフーリエ変換について


ルーチン等々はほぼ同じです。
下のプログラムは
\(
f(x,y)=\sin(4x)+\sin(3y)
\)
をフーリエ変換するものです。
フーリエ変換実行前の,位置空間での関数の形(“before.d”, x軸:\(x\), y軸:\(y\), z軸\(f(x,y)\)),
フーリエ変換し、波数空間で見た関数の形(“k.d”, x軸:\(k_x/(2\pi)\), y軸:\(k_y/(2\pi)\), z軸\(Imag\{g(k_x,k_y)\}\)),
逆フーリエ変換した後の関数の形(“after.d”, x軸:\(x\), y軸:\(y\), z軸\(f(x,y)\))は
before_k_after_2D_c
となります。
波数空間で\(k_x\)だけを見た時,\(k_x=0\)にあたかもピークが見えます。
これはy成分があるために出てくるもので、証明は以下のようにできます。
関数\(f(x,y)\)は
\(
\begin{align}
f(x,y)&=\sin(4x)+\sin(3y)\\
&=i\left\{-\frac{1}{2}e^{i4x}+\frac{1}{2}e^{-i4x}-\frac{1}{2}e^{i3y}+\frac{1}{2}e^{-i3y}\right\}
\end{align}
\)
なので、連続のフーリエ変換を考えるとすると、
\(
\begin{align}
&\int f(x,y)e^{-ik_x x}e^{-ik_y y}dxdy\\
&=i\left\{-\frac{1}{2}\int e^{i(4-k_x)x} dx \int e^{-ik_y y} dy+\frac{1}{2}\int e^{-i(4+k_x)x} dx \int e^{-ik_y y} dy \right. \\
&\;\;\;\;\;\;\;\;\;\;\;\; \left.-\frac{1}{2}\int e^{-ik_x x} dx \int e^{i(3-k_y)y} dy+\frac{1}{2}\int e^{-ik_x x} dx \int e^{-i(3+k_y)y} dy\right\} \\
&=i\pi\left\{-\delta(k_x-4)\delta(k_y)+\delta(k_x+4)\delta(k_y)-\delta(k_x)\delta(k_y-3)+\delta(k_x)\delta(k_y+3)\right\}
\end{align}
\)
となります。なので、フーリエ変換によってピークが出てくる位置というのは、
\((k_x,k_y)=(4,0), (-4,0), (0,3), (0,-3)\)
の4点となります。なので、y成分がいるために\(k_x=0\)にあたかもピークがあるように見えてしまうんですね。(証明終わり)

2次元のフーリエ変換はそのまま2次元データの配列を渡すことは出来ないため、一度、1次元配列への変換(サブルーチンswap1d2d)を呼び出しています。本来は2次元データと1次元データの間の関係ではEquivalenceを使うべきなのですが、何をしているのかの(個人的な)見やすさを考慮してわざわざ別の配列を用意して格納しています。

include "/opt/intel/mkl/include/mkl_dfti.f90"
program main
  implicit none
  integer::i,j,Nx,Ny
  double precision,allocatable::x(:),kx(:),mkx(:)
  double precision,allocatable::y(:),ky(:),mky(:)
  double precision::hx,xmin,xmax
  double precision::hy,ymin,ymax
  complex(kind(0d0)),allocatable::z(:,:),mz(:,:)
 
  double precision,parameter::pi=dacos(-1.d0)
  complex(kind(0d0))::func
  external::func
 
  Nx=100; Ny=100
  allocate(x(0:Nx-1),kx(0:Nx-1),mkx(0:Nx-1))
  x=0d0; kx=0d0; mkx=0d0
  allocate(y(0:Ny-1),ky(0:Ny-1),mky(0:Ny-1))
  y=0d0; ky=0d0; mky=0d0
  allocate(z(0:Nx-1,0:Ny-1),mz(0:Nx-1,0:Ny-1))
  z=dcmplx(0d0,0d0); mz=dcmplx(0d0,0d0)
 
  xmin=-3d0*pi
  xmax=3d0*pi
  hx=(xmax-xmin)/dble(Nx)
  ymin=-3d0*pi
  ymax=3d0*pi
  hy=(ymax-ymin)/dble(Ny)
 
  do i=0,Nx-1
     x(i)=xmin+dble(i)*hx
  enddo
  do j=0,Ny-1
     y(j)=ymin+dble(j)*hy
  enddo
  do i=0,Nx-1
     do j=0,Ny-1
        z(i,j)=func(x(i),y(j))  
     enddo
  enddo
 
  open(21,file="before.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')x(i),y(j),dble(z(i,j)),dimag(z(i,j))
     enddo
     write(21,*)
  enddo
  close(21)

  call dftf(Nx,kx,hx); call dftf(Ny,ky,hy) ! get frequency
 
  call dft2d(Nx,Ny,z,"forward") ! forward dft

  call dfts2d(Nx,Ny,kx,ky,z,mkx,mky,mz) !sort frequency
  open(21,file="k.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')mkx(i),mky(j),dble(mz(i,j)),dimag(mz(i,j))
     enddo
     write(21,*)
  enddo
  close(21)
 
  call dft2d(Nx,Ny,z,"backward") ! backward dft
  open(21,file="after.d")
  do i=0,Nx-1
     do j=0,Ny-1
        write(21,'(4e20.10e3)')x(i),y(j),dble(z(i,j)),dimag(z(i,j))
     enddo
     write(21,*)
  enddo
  close(21)
 
  stop
end program main

function func(x,y)
  implicit none
  double precision::x,y
  complex(kind(0d0))::func
 
  func=dcmplx(dsin(4d0*y)+dsin(3d0*x),0d0)
 
  return
end function func

がメインのプログラムで、以下が上で呼び出しているサブルーチンの中身です。

subroutine dftf(N,f,h)
  implicit none
  integer,intent(in)::N
  double precision,intent(in)::h
  double precision,intent(out)::f(0:N-1)

  integer::i
  double precision::mf(0:N-1)
 
  if(mod(N,2).eq.0)then
     do i=0,N-1
        mf(i)=(dble(i+1)-dble(N)*0.5d0)/(dble(N)*h)
     enddo

     do i=0,N-1
        if(i.le.N/2-2)then
           f(i+N/2+1)=mf(i)
        else
           f(i-N/2+1)=mf(i)
        endif
     enddo
  else
     do i=0,N-1
        mf(i)=(dble(i)-dble(N-1)*0.5d0)/(dble(N)*h)
     enddo

     do i=0,N-1
        if(i.le.(N-1)/2-1)then
           f(i+(N-1)/2+1)=mf(i)
        else
           f(i-(N-1)/2)=mf(i)
        endif
     enddo
  endif
 
  return
end subroutine dftf
!--------------------------------
subroutine dft2d(Nx,Ny,z,FB)
  !sikinote
  !date : 2015/12/29
  !developer : sikino
  !CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/deed.ja)
  use MKL_DFTI
  implicit none
  integer,intent(in)::Nx,Ny
  complex(kind(0d0)),intent(inout)::z(0:Nx-1,0:Ny-1)
  character(*),intent(in)::FB

  complex(kind(0d0))::z1(0:Nx*Ny-1)
  integer::Status
  TYPE(DFTI_DESCRIPTOR),POINTER::hand
  integer::id(1:2)
 
  !DFT : Discrete Fourier Transform
  !
  !n    --> number of data.
  !z(i,j) --> value of data at x and y
  !FB   --> "forward"  : Forward DFT
  !     --> "backward" : Backward DFT
     
  call swap1d2d(nx,ny,z,z1,1)
  id(1)=nx; id(2)=ny
 
  Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,2,id)
  Status = DftiCommitDescriptor(hand)

  if(trim(FB).eq."forward")then
     Status = DftiComputeForward(hand,z1)
  elseif(trim(FB).eq."backward")then
     Status = DftiComputeBackward(hand,z1)
  else
     write(6,*)"DFT string different"
     stop
  endif
  Status = DftiFreeDescriptor(hand)
 
  call swap1d2d(nx,ny,z,z1,-1)
 
  if(trim(FB).eq."backward")then
     z=z/dble(nx*ny)
  end if
   
  return
end subroutine dft2d
!---------------------------------------------
subroutine swap1d2d(nx,ny,z2,z1,isign)
  !change array matrix between 1D and 2D
  !  if isign ==  1  --> from 2D to 1D
  !  if isign == -1  --> from 1D to 2D
  implicit none
  integer,intent(in)::nx,ny,isign
  complex(kind(0d0)),dimension(0:nx-1,0:ny-1),intent(inout)::z2
  complex(kind(0d0)),dimension(0:nx*ny-1),intent(inout)::z1
  integer::j1,j2,k

  if(isign.eq.1)then
     do j2=0,ny-1
        do j1=0,nx-1
           k=j2*nx+j1
           z1(k)=z2(j1,j2)
        enddo
     enddo
  elseif(isign.eq.-1)then
     do j2=0,ny-1
        do j1=0,nx-1
           k=j2*nx+j1
           z2(j1,j2)=z1(k)
        enddo
     enddo
  endif
     
  return
end subroutine swap1d2d
!----------------------------
subroutine dfts2d(Nx,Ny,fx,fy,z,mfx,mfy,mz)
  implicit none
  integer,intent(in)::Nx,Ny
  double precision,intent(in)::fx(0:Nx-1),fy(0:Ny-1)
  complex(kind(0d0))::z(0:Nx-1,0:Ny-1)
  double precision,intent(out)::mfx(0:Nx-1),mfy(0:Ny-1)
  complex(kind(0d0)),intent(out)::mz(0:Nx-1,0:Ny-1)
  complex(kind(0d0))::mmz(0:Nx-1,0:Ny-1)
  integer::i,j,k,l  


  if(mod(Ny,2).eq.0)then
     do i=0,Ny-1
        if(i.le.Ny/2)then
           j=i+Ny/2-1
           mfy(j)=fy(i)
           mz(0:Nx-1,j)=z(0:Nx-1,i)
        else
           j=i-Ny/2-1
           mfy(j)=fy(i)
           mz(0:Nx-1,j)=z(0:Nx-1,i)
        endif
     enddo
  else
     do i=0,Ny-1
        if(i.le.(Ny-1)/2)then
           j=i+(Ny-1)/2
           mfy(j)=fy(i)
           mz(0:Nx-1,j)=z(0:Nx-1,i)
        else
           j=i-(Ny-1)/2-1
           mfy(j)=fy(i)
           mz(0:Nx-1,j)=z(0:Nx-1,i)
        endif
     enddo
  endif

  mmz=mz

  if(mod(Nx,2).eq.0)then
     do k=0,Nx-1
        if(k.le.Nx/2)then
           l=k+Nx/2-1
           mfx(l)=fx(k)
           mz(l,0:Ny-1)=mmz(k,0:Ny-1)
        else
           l=k-Nx/2-1
           mfx(l)=fx(k)
           mz(l,0:Ny-1)=mmz(k,0:Ny-1)
        endif
     enddo
  else
     do k=0,Nx-1
        if(k.le.(Nx-1)/2)then
           l=k+(Nx-1)/2
           mfx(l)=fx(k)
           mz(l,0:Ny-1)=mmz(k,0:Ny-1)
        else
           l=k-(Nx-1)/2-1
           mfx(l)=fx(k)
           mz(l,0:Ny-1)=mmz(k,0:Ny-1)
        endif
     enddo
  endif  

  return
end subroutine dfts2d

昔のMKLの離散フーリエ変換ルーチンについて


MKLは汎用性を高くするために仕様変更が行われ、
以前は

call zfft1d(...)

とするだけで順/逆方向フーリエ変換は使えていたようですが、いつからか(2004~2006年辺り?)廃止され、使えなくなっています。

間違い、勘違いがあると思います。
いかなる問題が発生しても私は責任を一切負いません。
それを念頭に置いたうえ使用してください。

参考文献


インテル®マス・カーネル・ライブラリーリファレンスマニュアル(2006年)
ニューメリカルレシピinC

スポンサーリンク

サバゲー用品まとめ

銃本体やパーツを安く、いろいろ見たい場合は
ガンショップFIRST
が良いと思います。

装備品の種類の豊富さ(服装、ポーチ等)では
S&Graf , エスアンドグラフ ミリタリー&カジュアルショップ
lro1
が群を抜いています。

人が少なくて、お店が綺麗で、見やすく、実際に銃に触れるお店は、
41PX 秋葉原
s41PXakihabara
がお勧めです。

購入したもの欲しいもの


銃本体

何もカスタムしなければ、0.5ジュール程度出ます。


mp7a1用マガジン

東京マルイ純正の多段マガジンは弾詰まりが頻発するそうです。使ってみた感想では、これは弾詰まりは起こりませんでした。


スコープ

初心者にはもってこいのものだと思います。
実際、工具なしでサイトの調節ができるのでとても使いやすいです。
前面にBB弾からスコープを守るカバーもついているので安心です。


mp7a1用バッテリー

【電動コンパクトマシンガン】対応バッテリーと充電器に関するご案内より、mp7a1用のバッテリー『No125 7.2V マイクロバッテリーEX CM-02』の生産は2015年5月に終了しています。
今後は在庫が尽きるまで提供は続けられますが、変換アダプター

を購入し、
No16 7.2V マイクロ500バッテリー
を購入することが推奨されています。




サバイバルゲームにはトレッキングシューズが最適かと思います。サバゲだけで履くのはもったいないので普段履きもできるこれがよさそうです。


無線

これは欲しいものです。やっぱりアンテナが長いほうが遠くまで聞こえるそうです。
ついでにこれも。


ヘッドセット
http://www.amazon.co.jp/gp/product/B001T7QJ9O?psc=1&redirect=true&ref_=ox_sc_act_title_1&smid=A2ZUZH56GAXKNY

無線があったらヘッドセットもそろえたいよね。

球体の空気抵抗と係数

まとめ


半径\(R\)の完全な球体に働く空気抵抗力\(F_d\)は
\(
\displaystyle F_d = 6\pi R \eta v + \frac{1}{2} C_d \rho \pi R^2 v^2
\)

です。方向も含めるのであれば
\(
\displaystyle \vec{F}_d = -6\pi R \eta |\vec{v}| \frac{\vec{v}}{|\vec{v}|} – \frac{1}{2} C_d \rho \pi R^2 |\vec{v}|^2 \frac{\vec{v}}{|\vec{v}|}
\)

となります。

ここで、

  • \(\eta\ \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) : 空気の粘性率(Viscosity) ※動粘度ではないです。
  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 空気の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 球体の速度の大きさ
  • \(R\ \mathrm{[m]}\) : 球の半径

を表します。
\(C_d\)は抗力係数(Drag Coefficient)と呼ばれ、厳密にはレイノルズ数\(R_e\)の関数\(C_d=C_d(R_e)\)となります。
レイノルズ数\(R_e\)とは無次元量で、
\(\displaystyle R_e=\frac{v \rho L}{\eta}\)
という量です。ここで\(L\)は、物体の大きさを表し、半径\(R\)の球の場合は直径\(L=2R\)に相当する量です。

球の場合の抗力係数\(C_d(R_e)\)をフィッティングしたものは、論文[2]より、
\(
\displaystyle \small C_d=\frac{24}{R_e}+\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
+ \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}} + \left(\frac{R_e^{0.80}}{461000}\right)
\)

と表されます。
Cdグラフ

20℃で1気圧(101325[Pa])、湿度0%の空気中では、
\(
\begin{align}
\eta &= 18.2\times 10^{-6}\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]} \\
\rho &= 1.205\mathrm{[kg \cdot m^{-3}]}
\end{align}
\)

という値になります。

温度、圧力、湿度を考慮する場合

温度\(T[\mathrm{K}]\)、圧力\(P[\mathrm{Pa}]\)の時,
空気の粘性率\(\eta\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\)は、
\(\displaystyle \eta\mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}=1.487 \times 10^{-6}\cdot \left(\frac{T\mathrm{[K]}^{1.5}}{T\mathrm{[K]}+117}\right)\)
であり、
空気の密度\(\rho\ \mathrm{[kg \cdot m^{-3}]}\)は、
\(\displaystyle \rho\ \mathrm{[kg \cdot m^{-3}]}=0.0034856447\cdot \frac{P\mathrm{[Pa]}}{T\mathrm{[K]}-0.670}\)
です。

この2つの式を[1]が正しいと思って確かめた温度の範囲-10℃~40℃で大体3桁くらい一致します。
この式は下に説明してある空気の場合を代入し、求めたものです。

2016/03/13 追)
湿度がある場合、空気の密度\(\rho_{w}\mathrm{[kg \cdot m^{-3}]}\)は湿度0%の時と比べて減少します。
温度\(t\mathrm{[{}^{\circ}C]}\)、湿度M(湿度60%の時、\(M=0.60\))、気圧\(P\mathrm{[Pa]}\)の場合、空気の密度\(\rho_{w}\mathrm{[kg \cdot m^{-3}]}\)は
\(
\displaystyle \rho_{w}\mathrm{[kg \cdot m^{-3}]}=
\left(0.0034856447\cdot\frac{P\mathrm{[Pa]}}{t\mathrm{[{}^{\circ}C]}+272.48}\right)\cdot
\left(1-M\cdot 0.378\cdot\frac{100\cdot 6.1078\times 10^{\frac{7.5\cdot t\mathrm{[{}^{\circ}C]}}{t\mathrm{[{}^{\circ}C]}+237.3}}}{P\mathrm{[Pa]}}\right)
\)

と記述されます。
ここで用いた関係式は
[5]に圧力pの水蒸気を含んだ空気の密度の関係式

[6]、[7]の飽和水蒸気圧を求めるTetensの式
を用いています。また、atmによる単位系では[8]に記述があります。


スポンサーリンク



スポンサーリンク

風の抵抗力の詳細


これと同じことは弾道計算(BB弾)の理論に記述してあります。
その中の風の抵抗力を部分を抜き出したものになります。

風(流体)の抵抗力の大きさ\(F_d\)は一般には次元解析により
\(
\displaystyle F_d=\frac{\eta^2}{\rho}\sum_{n=1}^{\infty}K_n \left(\frac{v \rho L}{\eta}\right)^n
\)

として与えられます[3]。ここで

  • \(\eta\ \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) : 流体の粘性率(Viscosity) ※動粘度ではないです。
  • \(\rho\ \mathrm{[kg \cdot m^{-3}]}\) : 流体の密度
  • \(v\ \mathrm{[m\cdot s^{-1}]}\) : 物体の速度の大きさ
  • \(L\ \mathrm{[m]}\) : 物体の大きさ(半径\(R\)の球の場合、\(L\)は直径\(L=2R\)に相当する量)
  • \(\frac{v \rho L}{\eta}\) : レイノルズ数(\(R_e\),Reynolds number)と呼ばれる無次元量

です。次元解析から具体的な定数の値\(K_n\)を求めることはできません。

速度の3乗以上に比例する項というのは見かけません。
おそらく、物体の動きを記述するために必要なパラメータがこの2つを考慮すれば含まれるため、\(n=3\)以上の項は冗長になるからでしょう。つまり、第一項には空気の密度\(\rho\)が含まれておらず、第二項まで入れることで\(\eta,\rho\)があらわに含まれるようになるからだと思います。

\(n=3\)以上の高次の効果は、抗力係数\(C_d\)に押し込めていると思います。

上の式の\(n=1,n=2\)をそのまま書き下すと
\(F_d=K_1\eta L v +K_2\rho L^2 v^2\)
となりますが、慣例として同じ次元を持つように式を変形させて
\(F_d=C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right)\)
の形で記述されます。ここでSは速度に垂直な、物体の断面積です。

\(C_1\)は物体が完全な球で、小さく、その速度がゆっくりである場合(レイノルズ数が2以下くらい)、流体の運動を記述する方程式、
ストークス方程式(Stokes’ law、ナビエ・ストークス方程式のレイノルズ数が小さい極限の方程式)
から導くことができて\(C_1=3\pi\)であることが理論的に導けます。

\(C_d\)は抗力係数(Drag Coefficient)と呼ばれ、
これはレイノルズ数の関数となり、
\(C_d=C_d(R_e)\)
となります。
この抗力係数は非常に厄介らしいです。物体の形状に依存したりするらしく、球の場合でも、どうやら実験的に求められているようです。
Wikipediaの抗力係数のページに若干の記述があります。
球の場合の抗力係数の関数\(C_d(R_e)\)をフィッティングした関数は論文[2]より、
\(
\displaystyle \small C_d=\frac{24}{R_e}+\frac{2.6\left(\frac{R_e}{5.0}\right)}{1+\left(\frac{R_e}{5.0}\right)^{1.52}}
+ \frac{0.411\left(\frac{R_e}{263000}\right)^{-7.94}}{1+\left(\frac{R_e}{263000}\right)^{-8.00}} + \left(\frac{R_e^{0.80}}{461000}\right)
\)

として書けるそうです。実際にプロットしてみるとこんな感じです。
Cdグラフ
確かに抗力係数のグラフと同じようになります。

よって半径Rの球の場合、風の抵抗力\(F_d\)の大きさは結局、
\(
\begin{align}
F_d &= C_1 L\eta v +C_d S \left(\frac{1}{2}\rho v^2\right) \\
&= 3\pi \cdot 2R\cdot \eta v +\frac{1}{2} C_d \rho \pi R^2 v^2 \\
&= 6\pi R \eta v + \frac{1}{2} C_d \rho \pi R^2 v^2
\end{align}
\)
として与えられます。

各定数の値
名称 捕捉
空気の粘性率(Viscosity)\(\eta\) \(\small 18.2\times 10^{-6} \\ \small \mathrm{[kg \cdot m^{-1}\cdot s^{-1}]}\) 空気中20度において。また、粘性率は数十気圧から数気圧の間圧力には依存しません。
が、温度に依存し、ある温度\(T_1[K]\)での粘性率\(\eta_1\)が分かっているならば、
温度\(T_2[K]\)での粘性率\(\eta_2\)はサザランドの定数Cを用いて、
\(
\eta_2=\eta_1 \left(\frac{T_1[K]+C}{T_2[K]+C}\right)\left(\frac{T_2[K]}{T_1[K]}\right)^{3/2}
\)
と記述されます。空気の場合、C=117です[4]。
乾燥した空気の密度\(\rho\) \(\small 1.205\mathrm{[kg \cdot m^{-3}]}\) 乾燥した空気の密度は温度と圧力に依存します。温度\(t[℃]\)、圧力をH[torr]とすると
\(
\rho=\frac{1.293}{1+0.00367t[℃]}\cdot \frac{H[torr]}{760}
\)
です[5]。

参考文献


[1]水・空気の物性 密度 粘度 動粘度
[2]Faith A. Morrison, “Data Correlation for Drag Coefficient for Sphere,” Department of Chemical Engineering, Michigan Technological University, Houghton, MI
[3]伊東敏雄著『な~るほど!の力学』学術図書出版社 (1994) p.220
[4]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.378
[5]国立天文台編『理科年表 平成21年 第82冊』丸善株式会社 (2009) p.376
[6]相対湿度の月別平年値 -理科年表オフィシャルサイト
[7]飽和水蒸気量 -wikipedia
[8] 空気 -wikipedia

ちなみに、
\(1 [\mathrm{torr}] = 133.322368 [\mathrm{Pa}]\\ 1[\mathrm{Pa}]=1[\mathrm{N\cdot m^{-2}}]=1[\mathrm{kg \cdot m^{-1}\cdot s^{-2}}]\)
という関係です。

fortran90でファイルの存在を確かめる

fortran90です。
ファイルが存在するかしないかを確かめるためのプログラムです。
組み込み関数”access”を用いると簡単に確認できます。

program main
  implicit none
  integer::access
 
  write(6,*)access( "./temp.d", " ")
 
  stop
end program main

上の文は、もしもファイル “./temp.d” が今いるディレクトリに存在するかどうかを返すものです。

また、access関数はファイルの有無(引数は空白文字” “)のみならず、そのファイルの読み取り権限(引数は”r”)、書き込み権限(引数は”w”)、実行権限(引数は”x”)を調べることができます。

ゼロが返ってきたらその権限(ファイルの有無)があるということです。
実行例は

 $ ls
main.f90  temp.d
 $ gfortran main.f90
 $ ./a.out
           0
 $

となります。intel®fortranコンパイラでももちろんできます。


スポンサーリンク


ここより以下のものはaccess関数を用いない方法です。


ファイルをopenできたら存在する、errorが返ってきたら存在しない、と判断します。
そのプログラムとサブルーチンはこんな感じ。

program main
  implicit none
  integer::ox
 
  call checkfile("temp.d",ox)
 
  write(6,*)ox
 
  stop
end program main

subroutine checkfile(fname,ox)
  implicit none
  character(*),intent(in)::fname
  integer,intent(out)::ox
 
  open(100,file=trim(fname),status='old',err=999)
  close(100)
  write(6,'(3A)')"file '",fname,"' exist"
  ox=1
  return
 
999 continue
  close(100)
  write(6,'(3A)')"file '",fname,"' don't exist"
  ox=0

  return
end subroutine checkfile

”ox”は〇×のつもりで使っています。
ファイル名”temp.d”をinputパラメータとして渡し、そのファイルが存在したらox=1を、存在しなかったらox=0を返します。

実行例は、

 $ ls
   temp.d main.f90
 $ gfortran main.f90
 $ ./a.out
   file 'temp.d' exist
           1
 $ gfortran main.f90
 $ ./a.out
   file 'temp1.d' don't exist
           0

となります。2回目の実行ではtemp.dをtemp1.dに変えて存在しない場合を表示しています。

スポンサーリンク

参考
ファイルの有無を調べる -知識の箱