第７回（2000年12月12日）講義録
〜高速フーリエ変換他〜

2次元配列

• 2次元配列 A(j,k)
• 例として行列の積を計算するプログラムを示す

       C   宣言部(寸法を10にしておく)
             integer n,j,k,m
             real A(10,10), B(10,10), AB(10,10)
       C   入力部
             open(10,file='matrix.data')
             read(10,*) n
             if ( n .gt. 10) then
                write(*,*) 'dimension is too large n =', n
                stop
             endif
             do  j = 1, n
                read(10,*) (A(j,k),k=1,n)
             enddo
             do  j = 1,n
                read(10,*) (B(j,k),k=1,n)
             enddo
       C   乗算部
             do  j = 1, n
  	      do k = 1, n
                   AB(j,k) = 0.0
                   do m = 1, n
                      AB(j,k) = AB(j,k) + A(j,m)*B(m,k)
                   enddo
                enddo
             enddo
       C   出力部
             do  j = 1, n
                write(*,*) (AB(j,k), k = 1,n)
             enddo
             stop
             end
       

DCLを使った高速フーリエ変換のプログラム

• 地球流体電脳倶楽部 DCL-5.0 を使って高速フーリエ変換を 行う。詳しくは配布のプリントを参照して。
• ここでは、簡易型フーリエ変換サブルーチンである、 EZFFTI, EZFFTF, EZFFTB を利用する。
• フーリエ変換の結果をmplotのスタイルで出力するプログラム。

        integer    N
        real       PI
        parameter (N=64, NN=N/2, PI=3.14159)
        integer    j
        real       WSAVE(3*N+15),R(N),A0,A(NN),B(NN),t(N),q
  
        write(*,*) 'Input parameter q !!'
        read(*,*) q
  
        do  j = 1, N
           t(j) = REAL(j-1)/N
           R(j) = sin(2.0*PI*q*t(j))   !ここで適当なデータをR(j)に与える
        enddo
  
        call EZFFTI(N,WSAVE)           !初期化
        call EZFFTF(N,R,A0,A,B,WSAVE)  !フーリエ変換
  C
  C   mplot用の出力
  C
        write(10,*) NN, 2              !データ数、グラフの本数
        do  j = 1, NN
           write(10,*) j, A(j), B(j)
        end do
        stop
        end
  

• 課題1: 次の関数をフーリエ変換しのグラフを表示せよ。ただし、データ数 N は128とせよ。

                       -1      (-1 < x < 0)
            f(x) =      0      ( x = 0)
                        1      (0 < x < 1)
  

• 課題2: 2次方程式の係数を読み込み解を求めるプログラムを作成せよ。