フォートランでの相互相関関数の計算

酒器について質問です。

ある２つのポイントで計測した実験データがあります。
その２つの実験データの間の時間のずれΔtを相互相関関数を用いて求めたいと考えています。

参考書やインターネットを通して理論は理解できるのですが、それをフォートランを用いた場合どのようにして計算させるのか分かりません。
具体的のどのようにして計算したらよいのでしょうか？
どのようなプログラム文を打てばよいか教えていただきたいです。

流れとしてはデータ数1000の場合
２つのファイルをopen
データ読み込み f()　g()
do p=-500,500
do i=1,1000
s=f(i)*g(i+p)
ss=ss+s
...........
のような流れになるのでしょうか？

どなたか分かる方是非教えてください
よろしくお願いします

No.1 ベストアンサー 回答者： irija_bari 回答日時： 2007/09/21 00:19

データ数が1000個ということは、

f(1), f(2), ..., f(1000)
g(1), g(2), ..., g(1000)
なのですよね？

ume3123さんのコードだと、g(i+p)の"i+p"のとる範囲が-499～1500になってしまいますよ。
相互相関関数をh(1000)とすると、

do p=式(p1), 式(p2)
h(p)=0.0
do i=式(i1), 式(i2)
h(p)=h(p)+f(i)*g(i+p)
enddo
enddo

とすればわかり易いのでは？
配列f(), g(), h()のインデックスが1～1000の範囲を逸脱しないように注意してください。

ヒント1：式(p?)には定数が入る。式(i?)にはpを使う（ただし片方だけ）。
ヒント2：1つのpに対して、h(p)=h(p)+f(i)*g(i+p)を1000回実行するわけではない。
(p=c(定数)の時、h(c)=h(c)+f(i)*g(i+c)を1000回実行するわかではないという意味です。)

がんばってください。

この回答へのお礼

ありがとうございます。
do p=1,10000
do i=1,10000-p
にしたらコンパイルできて一応出力でしました。

ただ本来ピークが真ん中にくるはずらしいのですが、端のほうにピークらしきものがあり，あかしな結果になりました。

本などを読んでいるとフーリエ変換しないとダメみたいなことが書いてあって、ますます分からなくなりました（＞＜）

実験でとった元のデータをdatファイルにして上記のようなプログラムでまわすだけは無理なんでしょうか？

お礼日時：2007/09/21 16:49

No.3 回答者： irija_bari 回答日時： 2007/09/21 20:22

ANo.1 です。

たびたびすみません。書き漏らしがありました。

>ただ本来ピークが真ん中にくるはずらしいのですが、

実験データを入れる前にプログラムが正しいことを確認しましたか？
あらかじめ時間のずれがわかっているデータ（例：sinとcos）で相互相関をとってみて、
プログラムが正しいことを確認すべきですよ。
実験はStep by Stepです。

No.2 回答者： irija_bari 回答日時： 2007/09/21 19:42

ANo.1 です。

==== Code ====
do p=式(p1), 式(p2)
h(p)=0.0
do i=式(i1), 式(i2)
h(p)=h(p)+f(i)*g(i+p-1)
enddo
h(p)=h(p)/(1000-p)
enddo
==== end of Code ====

Codeの下から2行目が抜けていました。ごめんなさい！！
この記述であってますかね？確認してください。
あと、g()もg(i+p-1)と記述すべきですね。ほぼ同じ結果になるはずですが。
indexが1から始まるプログラミング言語をしばらく扱っていない（0から始まるものが多い）
ため、他にも細かなミスがあるかもしれません。

do文については、ume3123さんの考え方でOKです。（細かなところまではわかりませんが。）
また、時間のずれを出すのに相互相関をとるという考え方も正しいはずです。

この回答へのお礼

ご丁寧にありがとうございます！
なんとかできました。

ただやはり波形処理を施さないといい結果は得れそうにないのでまた検討してみます！

お礼日時：2007/10/10 21:40