離散化フーリエ変換

　　（プログラムを参考にしてください。）

   0．Excel Macro   DFT   

最初にExcel Macroがあります。

Sheet1に

N,1,1,DT,

Time(i),X(i),Y(i)

のデータを貼り付け、実行すると、

Sheet2に計算結果が出力されます。

Sheet2には、周波数分析結果の

Power X(i),Fx(i),Power Y(i),Fy(i),TF(i),Phase(i)

が出力されます。

これをグラフにしたものを最初に図につけました。

macroのコピーのやり方。

新しいExcelのファイルを用意する。(***.xlsmで）

表示→マクロの表示→マクロ作成のフォームが

開くので、sub名をdftとして入れる。

マクロ作成画面になるので、End subを消して、

0.Excel macro   DFT  MACROと書かれた、

以下のマクロをコピーして（Sub dft()を除く）

貼り付ける。

Sheet1にデータを貼り付け実行すれば、

Sheet 2に計算結果が表示される。

Sheet1のデータ形式は、マクロを読んで

考えてね。

貼り付けるデータは、impuls_response.txtを

開いてSheet1の[A1]セルに貼り付ける。

カンマで区切られたデータです。

このデータの作成方法は、

FORTRANのimplus_response.forをコンパイル

して実行すると、impuls_response.txtができます。

試してみてね。固有振動数は1Hzです。

このマクロ（0．Excel macro  DFT  MACRO)は、

実験データを貼り付けて、実測データの周波数分析ができます。応用してね。

貼り付けるデータの1行目は、

N：データ数、1，1，dtです。

この、1，1、は意味はありません。

Nは、幾つにとっても、計算します。

FFTでは、ありません。DFTです。

Nを幾つにとってもいいところがDFTのいいところです。計算時間はFFTがNとするとDFTはN^2かかります。これが、難点ですが、最近のPCは、

計算速度が上がっているので、問題ありません。

時間軸と周波数軸で、それぞれ、切りのいい数値にしようとすると、ちょっと頭を使います。

謝辞

日野幹雄先生の本をだいぶ参考にさせて頂きました。良く書かれた本で、奥深いところまで書いてあります。ありがとうございました。

日野幹雄著　スペクトル解析　朝倉書店

YouTubeのヨビノリたくみさんの動画も参考に

させて頂きました。一部、動画で説明されていた

部分を清書させて頂き、使用しました。

フーリエ級数とフーリエ変換の関係がよくわかりました。ありがとうございました。

0．Excel macro        DFT  MACRO

      Sub dft()

      Dim x(16385) As Single      ' 時系列データ　x(t)

      Dim y(16385) As Single      ' 時系列データ　y(t)

      pai = 3.14159265358979

      Sheets("Sheet1").Select

      [A1].Select

      n = ActiveCell.Offset(0, 0)   '  n  データ数（時間領域のデータ数）

      dt = ActiveCell.Offset(0, 3)  '  dt 時間刻み

      For i = 1 To n

      x(i) = ActiveCell.Offset(i, 1)

      y(i) = ActiveCell.Offset(i, 2)

      Next i

      xave = 0#                    ' 初期化=0

      yave = 0#                    ' 初期化=0

      For i = 1 To n

      xave = xave + x(i)

      yave = yave + y(i)

      Next i

      xave = xave / n              ' x の直流分

      yave = yave / n              ' y の直流分

      For i = 1 To n

      x(i) = x(i) - xave           ' x の直流分カット

      y(i) = y(i) - yave           ' y の直流分カット

      Next i

      Tend = dt * (n - 1)          '  Tend : サンプリング時間の全体長

      f1 = 1# / Tend               '  f1 :  基本周波数　　:　　周波数分解能

      df = f1

      n2 = (n - 1) / 2             '  n2  :  周波数領域のデータ数

      dw = 2# * pai * df           '  dw :  基本角周波数

      w1 = dw

      Dim fxr(8192) As Single      ' x 　の周波数成分　の実数部

      Dim fxi(8192) As Single      ' x 　の周波数成分　の虚数部

      Dim fyr(8192) As Single      ' y 　の周波数成分　の実数部

      Dim fyi(8192) As Single      ' y 　の周波数成分　の虚数部

      Dim px(8192) As Single       ' x   のパワースペクトル

      Dim py(8192) As Single       ' y   のパワースペクトル

      Dim tfr(8192) As Single      ' 伝達関数の実数部   TF=y/x  x : input   y : output

      Dim tfi(8192) As Single      ' 伝達関数の虚数部   TF=y/x

      Dim tfa(8192) As Single      ' 伝達関数の振幅

      Dim phi(8192) As Single      ' 伝達関数の位相

      Dim QCOS(16384) As Single

      Dim QSIN(16384) As Single

      For m = 1 To n2              '    real part = ∫f(t)cos(mωt）dt    imaginary part = ∫f(t)(-sin(mωt))dt の m

      For i = 1 To n

      t = dt * (i - 1)             '    t  :  時刻

      QCOS(i) = Cos(m * w1 * t)    '    cos(mωt)

      QSIN(i) = -Sin(m * w1 * t)   '   -sin(mωt)

      Next i

      For mm = 1 To n2

      fxr(mm) = 0#                 '    初期化＝0

      fxi(mm) = 0#                 '    初期化＝0

      fyr(mm) = 0#                 '    初期化＝0

      fyi(mm) = 0#                 '    初期化＝0

      Next mm

      For i = 1 To n

      fxr(m) = fxr(m) + x(i) * QCOS(i)   ' ∫x(t)cos(mωt）dt

      fxi(m) = fxi(m) + x(i) * QSIN(i)   ' ∫x(t)(-sin(mωt))dt

      fyr(m) = fyr(m) + y(i) * QCOS(i)   ' ∫y(t)cos(mωt）dt

      fyi(m) = fyi(m) + y(i) * QSIN(i)   ' ∫y(t)(-sin(mωt))dt

      Next i

      fxr(m) = fxr(m) * dt

      fxi(m) = fxi(m) * dt

      fyr(m) = fyr(m) * dt

      fyi(m) = fyi(m) * dt

      px(m) = fxr(m) * fxr(m) + fxi(m) * fxi(m)

      py(m) = fyr(m) * fyr(m) + fyi(m) * fyi(m)

      tfr(m) = (fxr(m) * fyr(m) + fxi(m) * fyi(m)) / px(m)   ' 伝達関数の実数部

      tfi(m) = (fxr(m) * fyi(m) - fyr(m) * fxi(m)) / px(m)   ' 伝達関数の虚数部

      tfa(m) = Sqr(tfr(m) * tfr(m) + tfi(m) * tfi(m))        ' 伝達関数の振幅

      atfr = tfr(m)

      atfi = tfi(m)

      phi(m) = arctan(atfr, atfi)                            ' 伝達関数の位相

      Next m

      MsgBox (dt)

      Sheets("Sheet2").Select

      [A1].Select

      ActiveCell.Offset(0, 0) = "freq"

      ActiveCell.Offset(0, 1) = "power x"

      ActiveCell.Offset(0, 2) = "fx"

      ActiveCell.Offset(0, 3) = "power y"

      ActiveCell.Offset(0, 4) = "fy"

      ActiveCell.Offset(0, 5) = "tf"

      ActiveCell.Offset(0, 6) = "phi"

      For m = 1 To n2

      ActiveCell.Offset(m, 0) = df * m

      ActiveCell.Offset(m, 1) = px(m)

      ActiveCell.Offset(m, 2) = Sqr(px(m))

      ActiveCell.Offset(m, 3) = py(m)

      ActiveCell.Offset(m, 4) = Sqr(py(m))

      ActiveCell.Offset(m, 5) = tfa(m)

      ActiveCell.Offset(m, 6) = phi(m) * 180# / pai

      Next m

      MsgBox ("END")

      End Sub

      Function arctan(x, y)

      pai = 3.14159265358979

      If x = 0 And y = 0 Then

      arctan = 0

      ElseIf x > 0 And y = 0 Then

      arctan = ((pai / 2) * 0)

      ElseIf x = 0 And y > 0 Then

      arctan = ((pai / 2) * 1)

      ElseIf x < 0 And y = 0 Then

      arctan = ((pai / 2) * 2)

      ElseIf x = 0 And y < 0 Then

      arctan = (-(pai / 2))

      ElseIf x > 0 And y > 0 Then

      arctan = Atn(y / x)

      ElseIf x < 0 And y > 0 Then

      arctan = Atn(y / x) + pai

      ElseIf x < 0 And y < 0 Then

      arctan = -pai + Atn(Abs(y) / Abs(x))

      ElseIf x > 0 And y < 0 Then

      arctan = Atn(y / x)

      End If

      End Function

マクロは、ここまで。

以下は、FORTRANのプログラムです。

【説明のみ】

ここでの前提、MinGWをインストールしておくこと。

gfortran,C,C++をインストールのこと。

ネットで検索「MinGW」で出てくる。

ソース・リストは、このあと、下に記載。   

impuls_response.forは、1質点系のインパルスレスポンスを

求めるプログラムです。これが、DFTの入力データとなります。（例えばで用意しました。）

（任意のデータをDFTすれば、周波数分析が行えます）

04_DFT>gfortran impuls_response.for -o impuls_response.exe

   04_DFT>impuls_response

   04_DFT>gfortran DFT_16384.for -o DFT_16384.exe

  04_DFT>DFT_16384

C   実行の対話

C　ここでは、一自由度系のインパルス応答をフーリエ変換している。

C   

   04_DFT>DFT_16384

   INPUT FILE NAME

   impuls_response.txt

   Which Window

   1 : Rectangular Window

   2 : HANNING Window

   3 : HAMMING Window

   1

   OUTPUT FILE NAME OF SPECTRUM, TF, PHASE, COHR

   impuls_response-spectrum.txt

   OUTPUT FILE NAME OF AUTO-CORRELATION

   impuls_response-auto-correlation.txt

   OUTPUT FILE NAME OF CROSS-CORRELATION

   impuls_response-cross-correlation.txt

   04_DFT>

ここに出ているプログラムをコピー＆ペイスト

すると、1行おきにブランク行が出来るが、

構わず、******.forに保存して、コンパイル、

実行出来ることを確認しました。

FORTRAN　program List

　　(インパルス応答を求める、time,F(t),kx(t)を出力）

　（以下、Cの書いてある所から下をコピーしてTeraPadに

　　　貼り付け、impuls_response.forという名称で保存）

C   program file name  impuls_response.for

C

C   compile

C   >gfortran  impuls_response.for -o impuls_response.exe

C

C   実行

C   >impuls_response

C

C

C--------------------------------------------------------------

C

C   FORTRANのコンパイルには、MinGWを使用

C　　MinGWは、インターネットより入手可能。

C　　FORTRAN,C,C++をインストールのこと

C

C--------------------------------------------------------------

C

C       MAIN PROGRAM

      implicit real*8 (a-h,o-z)

      dimension f(3),q(3)

      common y(3),sm,sk,sc,t,force

      external sub1

      do 20 i=1,3

      y(i)=0.0

      q(i)=0.0

      f(i)=0.0

   20 continue

      ndt=8192

      dt=100./float(ndt-1)

      sm=10.0

      sk=394.7842

      zeta=0.01

      sc=zeta*2.0*sqrt(sk*sm)

      omega=sqrt(sk/sm)*(1.0-zeta*zeta)

      open(unit=3,

     &file='impuls_response.txt',

     &status='unknown')

      write(3,300)ndt,1,1,dt

  300 format(i5,2(',',i3),',',e16.8)

      t=0.0

      do 10 i=1,ndt

      t=float(i)*dt

      if(i.eq.1)then

      force=1.0

      else

      force=0.0

      endif

      call srkg(3,y,f,dt,q,sub1)

      write(3,100)t,force,y(2)*sk

  100 format(2(e16.8,','),e16.8)

   10 continue

      close(3)

      stop

      end

C

C

C     Runge Kutta Gill  methode

C

C     Variable

C     y(1)=t

C     y(2)=x

C     y(3)=dx/dt=v

C     

C     f(1)=dy(1)/dt=dt/dt=1.0

C     f(2)=dy(2)/dt=dx/dt=v=y(3)

C     f(3)=dy(3)/dt=d2x/dt2=dv/dt=α＝-c/m*dx/dt-k/m*x+force/m

C

C

      subroutine sub1(n1,f)

      implicit real*8 (a-h,o-z)

      dimension f(n1)

      common y(3),sm,sk,sc,t,force

      f(1)=1.0

      f(2)=y(3)

      f(3)=-sc/sm*y(3)-sk/sm*y(2)+force/sm

      return

      end

C

      subroutine srkg(n,y,ff,h,q,sub)

      implicit real*8 (a-h,o-z)

      dimension y(n),q(n),ff(n)

      dimension ccrkga(4),ccrkgb(4),ccrkgc(4)

      data ccrkga/0.5,0.292893218,1.70710678,0.166666666/

      data ccrkgb/1.0,0.292893218,1.70710678,0.333333333/

      data ccrkgc/0.5,0.292893218,1.70710678,0.5/

      do 10 i=1,4

      n1=n

      call sub(n1,ff)

      do 20 k=1,n

      xk=h*ff(k)

      qk=q(k)

      ry=ccrkga(i)*xk-ccrkgb(i)*qk

      s=y(k)

      y(k)=s+ry

      ry=y(k)-s

      q(k)=qk+3.0*ry-ccrkgc(i)*xk

   20 continue

   10 continue

      return

      end

FORTRAN　Program List（説明はコメント文にしているので大丈夫）

  (離散化フーリエ変換プログラム）

　（以下、Cの書いてある所から下をコピーしてTeraPadに

　　　貼り付け、DFT_16384.forという名称で保存）

C   program file name DFT_16384.for

C

C   compile

C   >gfortran DFT_16384.for -o DFT_16384.exe

C

C   exequte

C   >DFT_16384

C

C

C   実行の対話

C　ここでは、一自由度系のインパルス応答をフーリエ変換している。

C   

C   04_DFT>DFT_16384

C   INPUT FILE NAME

C   impuls_response.txt

C   Which Window

C   1 : Rectangular Window

C   2 : HANNING Window

C   3 : HAMMING Window

C   1

C   OUTPUT FILE NAME OF SPECTRUM, TF, PHASE, COHR

C   impuls_response-spectrum.txt

C   OUTPUT FILE NAME OF AUTO-CORRELATION

C   impuls_response-auto-correlation.txt

C   OUTPUT FILE NAME OF CROSS-CORRELATION

C   impuls_response-cross-correlation.txt

C

C   04_DFT>

C

C----------------------------------------------------------

C 注意書き

C  このプログラムを使用して製品を製作しても、計算結果を使用 

C　し、製品に不備が生じても当方は責任を負いません。

C　このプログラムは、無料公開するものですが、これで、商売を

C　行う事は、謹んで頂きたい。

C----------------------------------------------------------

C******************************************************

C*                                                    *

C*   SPECTRUM, COHERENCE AND CORRELATION BY DFT       *

C*   フーリエ級数に忠実にプログラミング　　　　　　　 *

C*                                                    *

C******************************************************

C

C     SPECTRUM IS OBTAINED AS ONE-SIDED SPECTRUM

C

C   X,Y         = GIVEN DATA

C   NMAX        = NUMBER OB DATA

C   DT          = TIME INTERVAL

C   NB          = POSITIVE POWER OF 2

C   NF          = NUMBER OF TERMS OF SMOOTHING FILTER IN FREQUENCY DOMAIN

C   PX1,PY1,PC1 = RAW SPECTRA OF X, Y AND CROSS-SPECTRUM  PC1

C   PX, PY, PC  = SMOOTHENED SPECTRA

C   PR1,  PI1   = REAL OR IMAGINARY PART OF CROSS-SPECTRUM  PC1

C   PR,   PI    = REAL OR IMAGINARY PART OF SMOOTHENED SPECTRUM PC

C   CXY,  CYX   = (OUT PUT) CROSS-CORELATION, I=+ OR I=-.

C   FX,   FY    = COMPLEX FOOURIER COMPONENT OF X,Y

C   FXR,  FXI   = COS OR SIN TRANSFORM OF X

C   FYR,  FYI   = COS OR SIN TRANSFORM OF Y

C   FRQ         = FREQUENCY

C

C------------------------------------------------------------------------

C

      CHARACTER FLIN*80,FLOUT*80,FLOUTAU*80,FLOUTCR*80,FLCHK*80

      COMMON /WINDOW/IWIND

      WRITE(6,100)

  100 FORMAT(1H ,'INPUT FILE NAME')

      READ(5,110)FLIN

  110 FORMAT(A)

      OPEN(UNIT=2,FILE=FLIN,STATUS='OLD')

   10 WRITE(6,*)'Which Window'

      WRITE(6,*)' 1 : Rectangular Window'

      WRITE(6,*)' 2 : HANNING Window'

      WRITE(6,*)' 3 : HAMMING Window'

      READ(5,*)IWIND

      IF(IWIND.LE.0.OR.IWIND.GE.4)GO TO 10

      WRITE(6,200)

  200 FORMAT(1H ,'OUTPUT FILE NAME OF SPECTRUM, TF, PHASE, COHR')

      READ(5,110)FLOUT

      OPEN(UNIT=3,FILE=FLOUT,STATUS='NEW')

      WRITE(6,*)'OUTPUT FILE NAME OF AUTO-CORRELATION'

      READ(5,110)FLOUTAU

      OPEN(UNIT=10,FILE=FLOUTAU,STATUS='NEW')

      WRITE(6,*)'OUTPUT FILE NAME OF CROSS-CORRELATION'

      READ(5,110)FLOUTCR

      OPEN(UNIT=11,FILE=FLOUTCR,STATUS='NEW')

C

C     入力データを読み込む（X(I)　、Y(I)　時系列データ）

C

      CALL READ5

C

C　　フーリエ変換をする。

C　　　　　生のスペクトルを求める。

C

      CALL FOURIE

C

C　　パワースペクトル、クロススペクトル、伝達関数、位相、コヒーレンス

C　　　を求める。

C

      CALL SPECTR

C

C　　　自己相関関数、相互相関関数を求める。

C

      CALL CORREL

      CLOSE(2)

      CLOSE(3)

      CLOSE(10)

      CLOSE(11)

      STOP

      END

C---------------------------------------------------------------------

C     入力データを読み込むサブルーチン

C　　　　　　　平均値（直流分）を求め、直流分をカットする。

C

C

      SUBROUTINE READ5

      COMMON /BL1/ X(16384),Y(16384)

      COMMON /BL4/ NMAX,DT,NB,NB2,NB4,NP,NF,NF2,MF,OMEG1

      READ(2,*)NMAX,NF,NF2,DT

      DO 10 I=1,NMAX

      READ(2,*)T,X(I),Y(I)

   10 CONTINUE

C      WRITE(3,600)

      SUMX=0.0

      SUMY=0.0

      DO 50 I=1,NMAX

      SUMX=SUMX+X(I)

      SUMY=SUMY+Y(I)

   50 CONTINUE

      XM=SUMX/FLOAT(NMAX)

      YM=SUMY/FLOAT(NMAX)

      DO 100 I=1,NMAX

      X(I)=X(I)-XM

      Y(I)=Y(I)-YM

  100 CONTINUE

C      WRITE(3,620)XM,YM

      RETURN

  500 FORMAT(3I5,G16.6)

  510 FORMAT(3G16.6)

  600 FORMAT(1H //10X,"I",12X,"X(I)",16X,"Y(I)")

  610 FORMAT(1H ,5X,I5,2E20.5)

  620 FORMAT(1H / 10X,"XMEAN=",E15.7,8X,"YMEAN=",E15.7// 10X,

     & "I",8X,"X(I)=X(I)-XM",8X,"Y(I)=Y(I)-YM"/)

      END

C--------------------------------------------------------------------

C      フーリエ変換をし、生のスペクトルを求める。

C

C　　　　　　　FXR(I) : Xの生のスペクトルの実数部

C              FXI(I) : Xの生のスペクトルの虚数部

C              FYR(I) : Yの生のスペクトルの実数部

C              FYI(I) : Yの生のスペクトルの虚数部

C

C

      SUBROUTINE FOURIE

      COMMON /BL1/ X(16384),Y(16384)

      COMMON /BL4/NMAX,DT,NB,NB2,NB4,NP,NF,NF2,MF,OMEG1

      COMMON /BL5/ FXR(8192),FXI(8192),FYR(8192),FYI(8192),S(4096)

      COMMON /WINDOW/IWIND

      COMPLEX XX(16384),YY(16384)

      DIMENSION QR(16384),QI(16384)

C

      NB=NMAX

C

C　　　ナイキスト周波数までのデータ（周波数）の数を求める。NB2

C

      NB2=NB/2

      NB4=NB/4

      PI=3.141592654

      TEND=DT*FLOAT(NMAX-1)

C

C　　　１周期（基本波）の時間の長さを求める。：TENDF

C

      TENDF=DT*FLOAT(NB-1)

C

C      WINDOW の設定　IWINDで場合わけ

C　　　１：Rectangular Window　　　（方形窓：何もせず）

C　　　２：HANNING Window

C　　　３：HAMMING Window

C

C

      WA=2.0*PI/TENDF

      IF(IWIND.EQ.2)GO TO 100

      IF(IWIND.EQ.3)GO TO 110

      GO TO 120

  100 DO 10 I=1,NB

      TX=FLOAT(I-1)*DT

      X(I)=X(I)*0.5*(1.-COS(WA*TX))

      Y(I)=Y(I)*0.5*(1.-COS(WA*TX))

   10 CONTINUE

      GO TO 120

  110 DO 20 I=1,NB

      TX=FLOAT(I-1)*DT

      X(I)=X(I)*(0.54-0.46*COS(WA*TX))

      Y(I)=Y(I)*(0.54-0.46*COS(WA*TX))

   20 CONTINUE

  120 CONTINUE

C

C      周波数分解能：DF

C

      DF=1./TEND

C

C　　　フーリエ変換の基本波の角周波数：OMEG1　　（ω）

C

      OMEG1=2.*PI*DF

C

C

C

      DO 1000 I=1,NB2

C

C      実数部の計算に使うQR(I)を作る。QR(I)＝COS（ωｔ）

C                                            COS(n ω1 t）

C　　　虚数部の計算に使うQI(I)を作る。QI(I)=－SIN（ωｔ）

C                                           -SIN(n ω1 t)

C

      DO 1020 J=1,NB

      T=DT*FLOAT(J-1)

      QR(J)=COS(OMEG1*FLOAT(I)*T)

      QI(J)=-1.*SIN(OMEG1*FLOAT(I)*T)

 1020 CONTINUE

C

      XXR=0.

      XXI=0.

      YYR=0.

      YYI=0.

      A=0.

C

C　　　フーリエ変換の式の計算。

C

C　　　　Xreal（ｆ）＝∫ｘ（ｔ）cos（ωｔ）ｄｔ

C

C　　　　Ximag（ｆ）＝∫ｘ（ｔ）（-sin（ωｔ））ｄｔ

C

C        X(J)とY(J)の両方について計算する。

C　　　　　（クロススペクトル、伝達関数、相互相関関数を求めるため）

C

C

      DO 1030 J=1,NB-1

      XXR=XXR+X(J)*QR(J)

      XXI=XXI+X(J)*QI(J)

      YYR=YYR+Y(J)*QR(J)

      YYI=YYI+Y(J)*QI(J)

 1030 CONTINUE

C

C　　　DTの掛算は最後に一括でする。

C

      FXR(I)=XXR*DT

      FXI(I)=XXI*DT

      FYR(I)=YYR*DT

      FYI(I)=YYI*DT

 1000 CONTINUE

      RETURN

      END

C-------------------------------------------------------------------

C    　パワースペクトル、クロススペクトル、位相、伝達関数、コヒーレンス

C　　　　　　PX(I)　　：　Xのパワースペクトル

C　　　　　　PY(I)　　：　Yのパワースペクトル

C　　　　　　PR(I)　　：　クロススペクトルの実数部

C　　　　　　PI(I)　　：　クロススペクトルの虚数部

C　　　　　　TFA(I)　 ：　伝達関数の振幅

C　　　　　　　　　　　　　TFA=Y/X　として計算

C　　　　　　　　　　　　　　　X　：　入力信号

C　　　　　　　　　　　　　　　Y　：　出力信号

C

C　　　　　　TFPHS(I) ：　位相(伝達関数、クロススペクトルとも同じ）

C　　　　　　COHR(I)  ：　コヒーレンス

C　　　　　　　　　　　　　入力のパワーがどれだけ出力側に伝わったかを

C　　　　　　　　　　　　　0.0～1.0の間の値で表現。

C

      SUBROUTINE SPECTR

      COMMON /BL3/ PX(8192),PY(8192),PR(8192),PI(8192)

      COMMON /BL4/ NMAX,DT,NB,NB2,NB4,NP,NF,NF2,MF,OMEG1

      COMMON /BL5/ FXR(8192),FXI(8192),FYR(8192),FYI(8192),S(4096)

      DIMENSION PX1(8192),PY1(8192)

      DIMENSION PR1(8192),PI1(8192),FRQ(8192),CPS(8192)

      DIMENSION PHASE(8192),COHR(8192)

      DIMENSION TFR(8192),TFI(8192),TFPHS(8192),TFA(8192)

C     *****  RAW SPECTRA AND CROSS- SPECTRUM  *******

      TN=2.0*DT*DT/(FLOAT(NB-1)*DT)

      OBT=1.0/(FLOAT(NB-1)*DT)*NF2

      TEND=DT*FLOAT(NMAX-1)

      RA=TEND/FLOAT(NMAX-1)

      BI=-TEND/FLOAT(NMAX-1)

      DF=1./TEND

C

C

C

      DO 50 I=1,NB2

C

C　　　　Pｘ＝（FXｒ ＋ i FXi）（FXr － i FXi）：共役複素数をかけて

C　　　　　　　　　　　　　　　　　　　　　　　パワースペクトルとする。

C

      PX(I)=(FXR(I)*FXR(I)+FXI(I)*FXI(I))

      PY(I)=(FYR(I)*FYR(I)+FYI(I)*FYI(I))

C

C       Pcross＝(FXr + i FXi)(FYr - i FYi)

C        上の式でクロスパワースペクトルを求める。

C　　　　下のプログラムでは、実数部と虚数部に分けて計算している。

C

C

      PR(I)=(FXR(I)*FYR(I)+FXI(I)*FYI(I))

      PI(I)=(-FXR(I)*FYI(I)+FXI(I)*FYR(I))

      CPS(I)=PR(I)*PR(I)+PI(I)*PI(I)

C

C       伝達関数：TFA　を求める。

C　　　　　TFA = (FYr + i FYi) / (FXr + i FXi)

C              下のプログラムでは、入力Xの共役複素数を分母、分子にかけて

C 　　　　　　 実数部：TFR、虚数部：TFIに分けて計算している。

C　　　　　　　この計算方法は、入力側のノイズに強い伝達関数が求まる。

C

C　　　　　　　出力側のノイズに強い伝達関数を求めるには、出力Yの

C　　　　　　　共役複素数を分母、分子にかけて求めればよい。

C

C

      BUNBO=FXR(I)*FXR(I)+FXI(I)*FXI(I)

      IF(BUNBO.EQ.0)THEN

      TFA(I)=0.0

      ELSE

      TFR(I)=(FYR(I)*FXR(I)+FYI(I)*FXI(I))/BUNBO

      TFI(I)=(FXR(I)*FYI(I)-FYR(I)*FXI(I))/BUNBO

      TFA(I)=SQRT(TFR(I)*TFR(I)+TFI(I)*TFI(I))

      ENDIF

C

C　　　位相：TFPHS　を求める。

C

      ATFR=TFR(I)

      ATFI=TFI(I)

      CALL SPHASE(ATFR,ATFI,BPHAS)

      TFPHS(I)=BPHAS

   50 CONTINUE

C     *****  COHERENCE AND PHASE     **********************

C　　　コヒーレンスを求める。

C

C

      PAI=3.141592654

      DO 100 I=1,NB2

      FRQ(I)=OBT*FLOAT(I)

      BUNBO=PX(I)*PY(I)

      IF(BUNBO.EQ.0.)GO TO 200

      COHR(I)=(PR(I)*PR(I)+PI(I)*PI(I))/(PX(I)*PY(I))

      GO TO 100

  200 COHR(I)=1.0

  100 CONTINUE

C     *****  OUT PUT  *************************************

C　　　計算結果の出力

C

C

      WRITE(3,625)

      WRITE(3,630)(I,FRQ(I),FXR(I),FXI(I),PX(I),FYR(I),FYI(I),

     &PY(I),PR(I),PI(I),CPS(I),TFA(I),

     &TFPHS(I)*180.0/PAI,COHR(I),I=1,NB2)

C     *****************************************************

      RETURN

  625 FORMAT("     I",

     &       "     FRQ(I)   ",

     &       "     FXR(I)   ",

     &       "     FXI(I)   ",

     &       "     PX(I)    ",

     &       "     FYR(I)   ",

     &       "     FYI(I)   ",

     &       "     PY(I)    ",

     &       "     PR(I)    ",

     &       "     PI(I)    ",

     &       "     CPS(I)   ",

     &       "     TF(I)    ",

     &       "     PHASE    ",

     &       "     COHR(I)  ")

  630 FORMAT(1H ,I5,13E14.6)

      END

C---------------------------------------------------------------------

C　　　位相：PHAを計算。

C

C

      SUBROUTINE SPHASE(BR,BI,PHA)

      PAI=3.141592654

      IF(BR.EQ.0.)GO TO 1200

      AB=BI/BR

      PHA=ATAN(AB)

      IF(BR.LT.0.AND.BI.GT.0.)PHA=PHA+PAI

      IF(BR.LT.0.AND.BI.LT.0.)PHA=PHA-PAI

      IF(BR.LT.0.AND.BI.EQ.0.)PHA=PAI

      GO TO 100

 1200 IF(BI.LT.0.)GO TO 1210

      PHA=PAI/2.

      IF(BI.EQ.0.)PHA=0.

      GO TO 100

 1210 PHA=-PAI/2.

  100 CONTINUE

      RETURN

      END

C--------------------------------------------------------------------

C　　　相関関数を計算。

C

      SUBROUTINE CORREL

      COMMON /BL3/ PX(8192),PY(8192),PCR(8192),PCI(8192)

      COMMON /BL4/ NMAX,DT,NB,NB2,NB4,NP,NF,NF2,MF,OMEG1

      DIMENSION TAU(4096),S(2048),TAU1(8192)

      DIMENSION CRX(4096),CRY(4096)

      DIMENSION CXY(4096),CYX(4096)

      DIMENSION W(8192),PW(8192),CRR(8192),CRI(8192)

C

C     ***** CORRELATION BY FFT OF SPECTRUM   *************

C

      NQ=NP-1

      NC=NB2

      NC2=NC/2

      NC4=NC/4

      DF2=1.0/(FLOAT(NB-1)*DT)

      TBDT=2.0*DT

C

C     逆フーリエ変換の計算

C

C　　　　自己相関関数は、パワースペクトルの逆フーリエ変換で求まる。

C　　　　　サブルーチンIDTFは、逆フーリエ変換の計算をする。

C

C　　　　相互相関関数は、クロススペクトルの逆フーリエ変換で求まる。

C　　　　　PCR　と　PR　は同じものである。COMMON文で同じ位置にある。

C　　　　　PCI　と　PI　も同じ。

C

      DO 110 I=1,8192

      PW(I)=0.0

  110 CONTINUE

      CALL IDFT(PX,PW,CRX,CRI,8192,4096,NC,NC2)

      DO 120 I=1,8192

      PW(I)=0.0

  120 CONTINUE

      CALL IDFT(PY,PW,CRY,CRI,8192,4096,NC,NC2)

      DO 130 I=1,8192

      PW(I)=0.0

  130 CONTINUE

      CALL IDFT(PCR,PW,CXY,CRI,8192,4096,NC,NC2)

      DO 140 I=1,8192

      PW(I)=0.0

  140 CONTINUE

      CALL IDFT(PW,PCI,CRR,CYX,8192,4096,NC,NC2)

C

      DO 100 I=1,NC2

      TAU(I)=TBDT*FLOAT(I-1)

      CRX(I)=CRX(I)*DF2

      CRY(I)=CRY(I)*DF2

      ST=(CXY(I)-CYX(I))*DF2

      CYX(I)=(CXY(I)+CYX(I))*DF2

      CXY(I)=ST

  100 CONTINUE

C

      DO 200 I=1,NC2

      IP=I+NC2-1

      IM=NC2-I+1

      CRR(IP)=CXY(I)

      CRR(IM)=CYX(I)

      TAU1(IM)=-TAU(I)

      TAU1(IP)=TAU(I)

  200 CONTINUE

      WRITE(10,700)

  700 FORMAT("     I",

     &       "     TAU      ",

     &       "     CRX      ",

     &       "     CRY      ")

      DO 10 I=1,NC2

      WRITE(10,710)I,TAU(I),CRX(I),CRY(I)

  710 FORMAT(I6,3E14.6)

   10 CONTINUE

      WRITE(11,720)

  720 FORMAT("     I",

     &       "     TAU      ",

     &       "  CROSS-CORREL")

      NCC2=NC2*2-1

      DO 20 I=1,NCC2

      WRITE(11,730)I,TAU1(I),CRR(I)

  730 FORMAT(I6,2E14.6)

   20 CONTINUE

      RETURN

C

  600 FORMAT(1H0//1H0,4X,"AUTO-CORRELATION  ,  CROSS-CORRELATION  OF

     &X - Y",/1H0,"  I ",5X,"TAU",8X,"CRX",8X,"CRY",7X,"CXY(+)",5X,

     &"CYX(-)")

  610 FORMAT(/(I4,1X,5E11.3))

C

      END

C-------------------------------------------------------------------

C  　　逆フーリエ変換

C

C　　　　Xｒ(ｔ）＝∫Fxr（ω）cos（ω t）df

C

C        Xi（ｔ）＝∫FXi（ω）sin（ω t）df

C

C                 Fxr 　: スペクトルの実数部　　　：　PX

C                 Fxi 　: スペクトルの虚数部　　　：　PW

C　　　　　　　　 Xr(t) : 時系列データ（相関関数）：　CRX

C　　　　　　　　 Xi(T) : 時系列データ（相関関数）：　CRI

C

C

      SUBROUTINE IDFT(PX,PW,CRX,CRI,ND,N2,N,NC2)

      COMMON /BL4/NMAX,DT,NB,NB2,NB4,NP,NF,NF2,MF,OMEG1

      DIMENSION PX(ND),PW(ND)

      DIMENSION CRX(N2),CRI(N2)

      DO 10 J=1,N2

      CRX(J)=0.0

      CRI(J)=0.

   10 CONTINUE

      DO 20 I=1,N

      DO 30 J=1,N2

      TT=2.0*DT*FLOAT(J-1)

      CRX(J)=CRX(J)+PX(I)*COS(OMEG1*FLOAT(I)*TT)

      CRI(J)=CRI(J)+PW(I)*SIN(OMEG1*FLOAT(I)*TT)

   30 CONTINUE

   20 CONTINUE

      PAI=3.14159265358979

      DO 40 J=1,N2

      CRX(J)=CRX(J)/2.0/PAI

      CRI(J)=CRI(J)/2.0/PAI

   40 CONTINUE

      END

【説明】

出力結果のspectrumの内容の入っているファイルを開いて、

DTRL+A、CTRL+Cですべてコピーし

次のExcelのマクロを入れたファイルのSheet1の[A1]セルに貼り付け、

実行し列を選択（周波数、TF(I)、・・)して散布図を描かせると

伝達関数のグラフが描ける。

ここでは、Sheet1に貼り付けたExcelのファイル名をマクロを有効とする、

xlsmにする。

そして、以下のマクロをマクロの関数に貼り付ける。

この伝達関数は、m=10,k=394.7842、zeta=0.01とした場合の

ラプラス変換で求めた伝達関数です。

[O1].selectで横に計算結果を並べています。

実行すれば、フーリエ変換とラプラス変換の比較が行えます。

マクロの下にその図を載せます。

Excel macro   sub transfer_function()

      Sub transfer_function()

      m = 10#

      k = 394.7842

      zeta = 0.01

      c = zeta * 2# * Sqr(m * k)

      omega0 = Sqr(k / m) * Sqr(1# - zeta * zeta)

      pai = 3.14159265358979

      f1 = omega0 / 2# / pai

      [B1].Select

      Dim freq(8193)

      For i = 1 To 4096

      freq(i) = ActiveCell.Offset(i, 0)

      Next i

      [O1].Select

      For i = 1 To 4096

      omega = freq(i) * 2# * pai

      bunbo = ((1# - (omega / omega0) ^ 2) ^ 2 + (2# * zeta * omega / omega0) ^ 2)

      tfr = (1# - (omega / omega0) ^ 2) / bunbo

      tfi = -(2# * zeta * omega / omega0) / bunbo

      tf = Sqr(tfr * tfr + tfi * tfi)

      theta = arctan(tfr, tfi)

      ActiveCell.Offset(i, 0) = freq(i)

      ActiveCell.Offset(i, 1) = tf

      ActiveCell.Offset(i, 2) = theta * 180# / 3.14159265358979

      Next i

      End Sub

      Function arctan(x, y)

      pai = 3.14159265358979

      If x = 0 And y = 0 Then

      arctan = 0

      ElseIf x > 0 And y = 0 Then

      arctan = ((pai / 2) * 0)

      ElseIf x = 0 And y > 0 Then

      arctan = ((pai / 2) * 1)

      ElseIf x < 0 And y = 0 Then

      arctan = ((pai / 2) * 2)

      ElseIf x = 0 And y < 0 Then

      arctan = (-(pai / 2))

      ElseIf x > 0 And y > 0 Then

      arctan = Atn(y / x)

      ElseIf x < 0 And y > 0 Then

      arctan = Atn(y / x) + pai

      ElseIf x < 0 And y < 0 Then

      arctan = -pai + Atn(Abs(y) / Abs(x))

      ElseIf x > 0 And y < 0 Then

      arctan = Atn(y / x)

      End If

      End Function