Fortranについて教えてください！

こんにちわ。
ヤマといいます。

以下について教えて頂けないでしょうか？

------------------- 途中省略 -------------------
DO 31 I = 1,NCARDS
LC = LC + 1
READ(8,FMAT) (XR(J), J=1,NPL)
IF(I.NE.NCARDS) GO TO 6
IF(JL.EQ.0) GO TO 6
JL=NPL+1-JL
DO 5 J=JL,NPL
5 XR(J)=0.
C---------------------------------------------------------------
C ONLY PRINT OUT A FEW OF THE FIRST AND THE LAST LINES of Input Motion
C---------------------------------------------------------------
6 ICHECK = NCARDS - I
IF (I .LE. 5 .OR. ICHECK .LT. 5) WRITE(6,2008) I,(XR(J), J=1,NPL)
IF (I .EQ. 10) WRITE (6,2009)
2009 FORMAT(3X,'........ INPUT MOTION READ NOT ECHOED...........')
C ENDIF
C
C FIND MAX. INPUT ACC. (XMAX)
C
311 DO 31 J = 1,NPL,2
N = N + 1
X(N) = CMPLX(XR(J),XR(J+1))
31 CONTINUE

...........................................

とあります。
いま、XRは、
----------------データ------------------------
-0.001694 -0.001668 -0.000086 -0.001356 -0.000678 0.000700 -0.001209 -0.000604
・・・・・・・・・・・・・・・・・・・・・
----------------データ------------------------
と1行に8列分のデータが存在します。
ソース中に「CMPLX(XR(J),XR(J+1))」とあります。
この意味がよく分かりません。通常だと実部にXR(J),虚部
にXR(J+1)を入れるというように解釈できるのですが、
本当でしょうか？ご存知の方ご意見よろしくお願いいたします。

No.9 回答者： semi-zzz 回答日時： 2005/08/09 14:18

＃１，４，６，７です。

ＦＦＴという情報が出ており、既に＃５，８で解答が出ているようです必要ないと思いますが、
もし、ＦＦＴのプログラミングを勉強中でしたら、以下の書籍にも理論及びFORTRANのソースが書かれていますので、参考にされてはいかがでしょうか？

日野　幹雄 著「統計ライブラリー スペクトル解析」1977年 朝倉書店　ISBN 4-254-12511-9

大崎　順彦 著 「新・地震動のスペクトル解析入門」1994年 鹿島出版会　ISBN 4-306-03270-1

後者はスペクトルの定義が一般のものと若干異なった定義をしているので取扱は注意する必要がありますが、フーリエ解析の入門書としてはお勧めです（建築・土木系の方でしたら特に）。

No.8 回答者： tatsumi01 回答日時： 2005/08/09 11:38

No. 2 = No. 5 ですが蛇足です。

実データ x(0), x(1), ...., x(N-1) の DFT は
F(k) = Σ[j=0, N-1] W^jk x(j), ここで W = exp(-i 2π/N): 回転因子, i は虚数単位
x(j) を偶数番目と奇数番目にわけると
F(k) = Σ[j=0, N/2-1] W^2jk x(2j) + Σ[j=0, N/2-1] W^(2j+1)k x(2j+1)
= Σ[j=0, N/2-1] U^jk x(2j) + W^k Σ[j=0, N/2-1] U^jk x(2j+1), ここで U = W^2: N/2 での回転因子
となって、最右辺第１項は偶数番目のデータの DFT、第２項は奇数番目のデータの DFT です。ただし、サンプル数は N/2 と半分になっています。
お尋ねの問題では、偶数番目 x(2j) を実部、奇数番目 x(2j+1) を虚部に入れて DFT (N/2) を計算するのですが、得られた DFT とその共役複素数から、実部、虚部それぞれの DFT (N/2) を簡単に計算できます。それを上の式に代入することにより、サンプル数 N の DFT が計算できます。N が大きくなると計算量は FFT の計算量が支配的で、ほぼ半分の計算時間ですみます。
昔の FORTRAN は配列の添え字に 0 が許されないので、DFT の式との対応で苦労したものです。

No.7 回答者： semi-zzz 回答日時： 2005/08/08 18:32

＃６の訂正

FORTRANの場合、累乗を示すのは
「＾」ではなくで「**」の間違いでした。

時刻歴波形を用いて複素数を計算するのはＦＦＴがぱっと思いつきます。
たしか1960年代半ばにＦＦＴを用いた計算法が用いられるようになったので、プログラムの作成時代からいっても、＃５さんの解答にもあるようにＦＦＴ（高速フーリエ変換）でしょうか？

この回答へのお礼

ご回答ありがとうございます。その通りです。FFTの計算に使用いたします。

お礼日時：2005/08/09 10:21

No.6 回答者： semi-zzz 回答日時： 2005/08/08 18:19

＃１，４です。

＃４を書いている間に＃２さんに補足していたようで、＃２さんへの補足で入力データが時刻歴データであることがわかりました。
入力データは加速度時刻歴ですね（質問文のコメントにACCとあるので）。

＃４で書いた内容と相反するのですが、なぜ複素数化するのかと言うことと、複素数の絶対値をとると言うことの意味を考えてみましたので、書き込みします。

実数上で扱うと、サブルーチンの計算部分は以下のような内容になります。
XMAX=0.
DO 1 I=1,MX
XA = SQRT(XR(I)^2+XR(I+1)^2)
IF (XMAX.GT.XA) GO TO 1
NXMAX = I
XMAX = XA
1 CONTINUE

これは連続した２つのデータの２乗和の平方根を求めています。
なぜそんなすること、計算目的が書いていないのではよくわかりませんが、エネルギーとしてのデータに変換する（速度データでしたら速度の２乗は運動エネルギーに比例する量なのですが、加速度のようですから物理的な意味はちょっとわかりません）、平均化・スムージングするなどの理由を思いつきます。
なお、もしXA = SQRT((XR(I)^2+XR(I+1)^2)/2)なら時間刻み分の時間を時定数とした実行値を求めていることになります。

わざわざ複素数化して計算したのは、計算上のテクニックだと推定されます。
特に補足にあったサブルーチンのコードされた年が1971年10月になっており、データ入力もカードから読み込んでいるような時代のもののような記述になっています(NCARDSとかいう変数があるので）。

このプログラムが書かれた当時のコンピュータは汎用の大型計算機でもかなりメモリが小さかったために、大量データの計算を行う際にメモリを有効に使う方法として一見変な方法を用いて計算するという手法を高度なプログラマーは用いていることがあったのでそのためではないかと思います。

ちなみにX(J)=CMPLX(XR(J),0.0)という複素変換は実数の時刻歴波形をから周波数分析を行う際に、複素領域上の計算が必要なFFTのプログラムなどでよく使用します。

都合により本日はこれにて失礼しますが、補足が必要でしたら書き込んでおいてください。明日なら補足できます。

No.5 回答者： tatsumi01 回答日時： 2005/08/08 18:14

No. 2 です。

私のアドバイスへのお礼で、「時時刻暦波形（意味不明）」を左から順に X に入れるとあります。
しかし、プログラムを見ると波形（恐らく等間隔にサンプリングされた値）が配列 XR に読み込まれ、それを２個づつまとめて複素数配列 X に格納しています。
これから推測されることはＦＦＴをさらに高速（２倍速）にするためではないかと思われます。普通にＦＦＴを行うには、実部にデータを入れ虚部を０としてデータ配列を作り、ＦＦＴサブルーチンを呼びます。しかし、これではデータが実数であるという特性を生かしていません。
データが実であるとき、この例のように実部と虚部に交互にデータを入れ、ＦＦＴを呼ぶとデータ数が半分で済みます（ＦＦＴの出力に細工をする必要はあります）。
それ以外に、同質のデータを実部と虚部に交互に入れる必要性は思い浮かびません。

この回答へのお礼

おはようございます！失礼いたしました。「時時刻暦波形」ではなく「時刻暦波形」の間違いでした。これは具体的に言えば、地震波形のデータを並べたものになります。
それと、複素数配列Xの実部と虚部にに実であるデータを交互に入れるという方法ならばデータ数が半分で済むということは初めて知りました。賢いですね、ありがとうございました。

お礼日時：2005/08/09 10:20

No.4 回答者： semi-zzz 回答日時： 2005/08/08 16:57

＃１です。

補足に対してアドバイスします。

計算目的がよくわからないので、推定で書きます。データ自体は時刻歴のようなもののデータのようですので、何の計算をしているのかもわかりません。
何のための計算をするプログラムなのかを書いてもらえればもう少し推測できるのですが、

まず、複素数の絶対値を求める組み込み関数ははABSではなく、CABSです。

つぎに、サブルーチンでdimension文でX(1)と型宣言していますので、サブルーチン内ではＸは実数となっています。またABSを用いていること、XMAX = 0.となっていることもサブルーチン内では、Ｘは実数扱いとなっています。すなわちサブルーチン内では単純に実数Ｘの絶対値を求めています。
データが時刻歴のような時間変動するデータでしたらこのサブルーチンはその時間内のデータ最大値を求めているだけです。

なぜこんなことになっているのかは、call文がきちんと書いてないので推定ですが、上記のことから、主文とサブルーチン内で同じ名称の変数Ｘを用いていても、同じ変数を意味していない可能性があります。すなわち変数名は一緒でも違う数値を取り扱っているということです。

call分の引数の位置とsubroutine文の引数の位置をチェックしてみてください。call文の最初の位置に来ている変数の値ががサブルーチン内のＸに引き継がれています。
SUBROUTINE XMX(X,MX,XMAX,NXMAX)

万が一主文で複素数のＸが実数としてサブルーチン内に引き継がれているようでしたら、計算自体うまくいかないものと思います。

なお、FORTRANではサブルーチン内と主文などで変数名は独立しており、同じ名称をしても関係のない変数として取り扱えます。主文とサブルーチン文の変数を関係付けているのはcall文とsubroutine文の並びやcommon文です（本ケースでは使用していないようですが）。

No.3 回答者： cubics 回答日時： 2005/08/08 16:17

なぜ複素数の絶対値を取るか、そもそも問題が書いてないので、決定的なことは何もわかりませんが、推測できるのは、データが２個ずつ意味をもつこと、たとえば平面座標系の座標値として、減点からの距離が最大な座標を求める問題、といったことですね。

複素数を用いたのは、組み込み関数で、計算のステップを減らしたということしか、これだけの質問では推測しかねます。

No.2 回答者： tatsumi01 回答日時： 2005/08/08 16:10

No. 1 への補足についてですが

「そもそもどうして、こんなことをするのか分かりません。」は、このサブルーチンの内容でしょうか。複素数の配列の中での絶対値の最大値と、その要素の番号を求めるサブルーチンですね。なお、X には Complex 宣言が必要です。「X(J)=CMPLX(XR(J),0.0)とすべきか」と書かれていますが、それでは複素数にはなりません。読み込んだデータを実部に入れただけで、本質的には複素数になっていません。このプログラムでは、入力データを２個１組にして、実部・虚部にかわるがわる入れた複素数を作っています。その複素数配列の中で絶対値の最大値が必要なんです。
もしかして、最初の質問のプログラムで
JL=NPL+1-JL
DO 5 J=JL,NPL
5 XR(J)=0.
としている部分のことでしょうか。"JL=NPL+1-JL" の右辺の JL は最初どこかで値が定められており、左辺の JL は後ろから数えた JL 番目の番号です。この DO ループは JL より後ろの値を全部０にクリアする目的です。（ただし、このように一つの変数に二つの意味を持たせるのはバグ要因です。）なぜ０でクリアする必要があるかというと、恐らく配列 X は繰り返し使われるので、データ数が少ないとき配列の後ろに以前の値が残ると困るので掃除しています。
なお、プログラムにはいくつかバグがあるように見受けられます。

この回答へのお礼

大変申し訳ありません。実はこのXR(J)に入る値は時時刻暦波形を左から右に順番に時間刻み0.02秒毎に示しています。つまり、X(J=0)は時刻0sの値、X(J=1)は時刻0.02秒での値、X(J=2)は時刻0.04秒での値、・・・・と与えています。したがってXRの値そのものが実数値になっています。ご指摘のように実部にXRを入れて、虚部に0を入れてもXは複素数にはなりません。ここでは文が長くなると投稿できないので、前文・後文を省略させて頂いておりますが、確かにXはComplexで宣言しています。JLの値も実は前文に出てくるので、ちゃんと入っています(JL=4)。

お礼日時：2005/08/08 16:51

No.1 回答者： semi-zzz 回答日時： 2005/08/08 15:27

CMPLXはFORTRANの組み込み関数で、複素数変換を行うためのものです。

パラメータを２つ取り、CMPLX(a,b)としたとき、質問文にあるとおり、aが実部、bが虚数部の数字になります。

また上記の例でいうと、最初に型宣言文COMPLEXでＸが複素数であることを定義しておく必要があります(FORTRANでは複素数は１つの変数として取り扱えます)。

ちなみに複素数の変数から実部を取り出すにはREAL,虚部を取り出すのにはAIMAGという組み込み関数を使用します。

この回答へのお礼

ご回答ありがとうございます。やはりそうですか。。。このXR(J)に入る値としましては、時間間隔が0.02秒毎に与えられている8F10.6形式の数字を読み込んでいきます。複素数X(N)に入れたものを今度は、
サブルーチンXMX(X,MA,XM,NXMAX)をコールして、それらXの内の最大値を取ろうとしています(入力:X,MA,出力:XM,NXMAX)。ちなみに、
C*************************************
SUBROUTINE XMX(X,MX,XMAX,NXMAX)
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C THIS ROUTINE FIND MAX. VALUE, XMAX, AND NUMBER OF MAX. VALUE, NXMAX.
C OF ARRAY X WITH MX NUMBER OF VALUES
C
C CODED PER B SCHNABEL OCT. 1971
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
DIMENSION X(1)
XMAX = 0.
DO 1 I = 1,MX
XA = ABS(X(I))
IF (XMAX.GT.XA) GO TO 1
NXMAX = I
XMAX = XA
1 CONTINUE
RETURN
END
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
です。
そもそもどうして、こんなことをするのか分かりません。
それをするならば、X(J)=CMPLX(XR(J),0.0)とすべきかかと思いました。なぜとなりの列を虚部にする必要があるのでしょうか？

お礼日時：2005/08/08 15:46