コヒーレンス関数について

質問

二つの波形の一致率を調べるための方法を聞いたところ「コヒーレンス」を使えばいいとのことで、今コヒーレンス関数を調べております。
http://oshiete1.goo.ne.jp/kotaeru.php3?q=2314537

で、一応式は検索できたのですが、
１、http://dspace.wul.waseda.ac.jp/dspace/bitstream/2065/396/9/Honbun-31_chapter3b.pdf
２、http://www.hulinks.co.jp/support/flexpro/v7/dataanalysis_crsp.html#02
３、http://www.cybernet.co.jp/matlab/support/manual/r14/toolbox/signal/spectra7.shtml
４、http://www.onosokki.co.jp/HP-WK/c_support/tech_term/cf_fft/cf3_2.htm
たぶん２，３，４は同じ事を書いてるのだと思うのですが、１は左辺に二乗があり、右辺に二乗がありません。これは

A)クロスパワースペクトル密度とクロススペクトル、クロスパワースペクトルは別物と言う事でしょうか？
B)ｒ＾２がクロススペクトルなのでしょうか？それともｒがクロススペクトルなのでしょうか？

次に、コヒーレンス関数の式が上の１，２，３，４で求められるとしたら、パワースペクトル、クロスパワースペクトルとはいったいどうやって計算すれば良いのでしょうか？

C)パワースペクトル・・・FFTの計算結果の√（実数部＾２+虚数部＾２）と言う事でいいんでしょうか？

D)クロススペクトル・・・FFTの計算結果で、XとYで（実数部、虚数部）がそれぞれ(X,iX)、(Y,iY)があった時、クロススペクトルは
（（X+iX）*（Y-Yi））＾２＝（XY-iXY+iXY+XY)^2 = (2XY)^2
で良いんでしょうか？（虚数部の計算は打ち消しあうので実数部だけでいいんですか？）

よろしくお願いします。

回答No.4

noname#65504 · ベストアンサー

＃１～３です。

＞データが2048個です。
＞もしかしたら1024か2048辺りで割った方がいいのでしょうか？

というと、Ｎ／２がかかっていそうですね（Ｎ：データ個数）。
このＮ＊Δｔ／２（Δｔ：サンプリング間隔）＝サンプリング時間／２という係数をかけてフーリエスペクトルとしているのは、東京大学の大崎先生が有名です（建築学）。
そして大崎先生は以下の著書で、フーリエ振幅スペクトルをそう定義して、fortranで書いたソースを公開しています。質問者が参考にしているのは、その流れをくむ計算アルゴリズムのように思います。

「地震動のスペクトル解析入門」　大崎　順彦著、鹿島出版会　

現在は改訂版が出て、「新・地震動のスペクトル解析入門」という本になっています（ISBN:4306032701）。でも旧版の方が説明が丁寧で良い本なのですが。

上記の本は地震工学の分野にしぼってフーリエ解析までしか説明していませんので、おそらく分野が異なる質問者が購入してもあまり意味がなさそうですが。もし建築・土木分野の人でスペクトル解析を勉強したいのでしたらお薦めの本です。

また、パワースペクトルを周波数で基準化したものをパワースペクトル密度と呼びます（「ディジタル信号処理入門」　木戸健一著　丸善　p20：小野測器はこの著書を参考文献として紹介しているので、小野測器とは用語の定義は同じと思います。ただし絶版のため入手困難です）

Ｎ＊Δｔ／２＝（Ｎ／２）＊（１／サンプリング周波数）ですから、パワースペクトルではなく、大崎先生の定義する（一般的でない定義）フーリエスペクトルから、パワースペクトル密度を計算しているようにも思います。
このあたりの影響がどうなるかはちょっとわかりません（式をきちんと追えばわかると思いますが）。

＞「スペクトル解析」日野幹雄著　の本を買ってきたんですが、それでも「？？？」でした。

日野先生の「スペクトル解析」はこの分野の定番の本ですが、かなり難しいですね。私もよくわからないことが多いです。初学者には敷居が高いので仕方ないと思います。

＞そして、何個の平均を取れば良いのでしょうか？

私がアナライザーを用いて解析する際には少なくとも最低でも１０個、通常１００個ぐらいのデータを使用しています。でも正弦波だとばらつきがないので、平均化しても結果が変わらないので、コヒーレンスの検証データとしては利用できないかもしれません。

＞この平均なんですが、どれの平均ですか？？　

平均化処理というのは、時間領域で平均化する方法もあるのですが、コヒーレンスを求める際にする平均化処理は、N個の信号からなるデータをM組用意して、ＦＴＴ変換したスペクトル上で、各周波数ごとに平均値を求めるものです（申し訳ありませんが、パワー、フーリエ、クロスのどの段階ですればよいかはよく覚えていません）

小野測器のホームページで見るのならＦＦＴアナライザーの使用法について解説してあるＰＤＦファイルをご覧下さい（ＰＤＦにリンク貼るのはこのサイトでは違反らしいので、貼れません）。

＞Ar[i] = (Ar[i] + Ar[i+1] + ar[i+2]) / 3

ちなみに上記の処理は平均化処理ではなく、複雑なスペクトルから全体的な傾向などを見るために行うスペクトルの「平滑化処理」です（平滑化処理の詳細は大崎先生の著書に書いてあります）。

お礼

大崎先生の流れかどうかは分かりませんが、参考にしたサイトはこちらです。
http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/
そして、残念ながら建築・土木分野の仕事ではないです。

平均化と平滑化は同じものだと思っていたのでとても参考になりました。

>ＦＴＴ変換したスペクトル上で、各周波数ごとに平均値を求めるものです
と言う事は数回測定したＦＦＴの結果を平均すると言う事ですよね？どうも以前ＤＬしたソフトのクロススペクトルは平滑化しているようでした。
http://www.ritsumei.ac.jp/se/rv/izuno/soft.html　コヒーレンスの計算

ただ、一度ＦＦＴしてからパワー・クロススペクトルを出すのと、相互・自己相関関数の後にＦＦＴするのだと後者の方が平滑化したようなデータになったので出来るかと思ったのですが、どうやっても無理でした。（１を超えるデータが出来てしまう事がある）

最終的にＤＳＰで信号処理を考えていたのですが、どうやら処理時間が間に合いそうも無いためコヒーレンスはあきらめようと思います。（ＦＦＴ１回でも処理時間が厳しいです）

長々と付き合ってもらった上に結果が出なくて申し訳ありません。
ありがとうございました。

回答No.3

noname#65504

＃１です。

>音・振動のスペクトル解析 　工博 金井浩　ISBN：4-339-01105-3

この本はかなり応用編ですね。
コヒーレンスは信号処理の中でも高度な部類ですので、コヒーレンスについて書いていなくとも信号処理の入門書を読んでから、コヒーレンスまで書いてある本に進んだ方がよいかもしれません。

また信号処理の入門書の後に出ている本はディジタル信号処理関係に比べて出版されている本がぐっと少なくなるのですが、モード解析の本を読んでみるのもよいかもしれません。コヒーレンスや伝達関数を使った解析法について述べてあることが多いです。

＞自分でやる場合ですが、ＦＦＴを求めるとＳＩＮ波なのである周波数が1024で、他は０になります。
次に、パワースペクトルですが、ある周波数が1024で他は全部０でした。

１０２４というのは、振幅（縦軸の値）ですか？
フーリエ振幅スペクトルの振幅と正弦波の振幅は関係があります。
著者によってフーリエ上の振幅に係数がかかっていることもあるのですが、基本的に正弦波の振幅とフーリエ振幅の値はほぼ同じになるはずです（計算誤差などがあるので一致はしないこともある）。
また振幅が最も大きくなる周波数と入力正弦波の振動数もほぼ一致するはずです。
ここはよくチェックしてみてください。

ちなみにＦＦＴでは解析に用いるデータ個数は２の累乗個のデータであることが必要なので、それが１０２４ではないでしょうか？
もし入力に用いた正弦波の振幅値が１なら振幅値にデータ個数がかかっているように思います。

また上記の状況から窓関数は使用していないのではないかと思います（窓関数がかかると正弦波の振動数付近にも０以外の値が現れるので）。
切り取った信号のがちょうどきりのいいところできれていれば正弦波なら問題ないのですが、複雑な信号を処理する際には目的に合わせた窓関数（ウィンドウ）をかけて処理しないと、リンク効果により変な結果が出ますので、窓関数を使用した方がよいです。
窓関数については小野測器のホームページにいろいろ解説がありますのでよく見てください。

＞小野測器のＨＰでは平均しなければ駄目と書いてましたし

コヒーレンスを求めるには平均化処理をしなければなりません。

ここは訂正です。
＞C)パワースペクトル・・・FFTの計算結果の√（実数部＾２+虚数部＾２）と言う事でいいんでしょうか

この分野では人によって多少言葉の定義が人によって異なるのですが、この式で表されるものは、私はフーリエ振幅スペクトルまたは単純にフーリエスペクトルと呼んでいます。そして平方根を取らないもの（２乗和）をパワースペクトルと呼んでいます。

＞他にも何かする事があるのかな？と思ってます。

平均化と窓関数のほか、入力信号をディジタル化する前に、アンチエリアジングフィルターというアナログフィルターをかける必要があります。

＞あと少しのところまで来てるのではないか？と思っているのですが今悩んでます。

お礼を見るとそんなところまで来ていそうですね。もう一息ですので頑張ってください。

お礼

ありがとうございます。
とりあえずテストで試してたのが振幅１のＳＩＮ波で、データが2048個です。（取り込んでるのではなく、計算で出したＳＩＮ波です）
窓関数を使うと窓関数の周波数も入ってしまい、比較しにくくなるので使ってませんでした。

ＦＦＴは以前もやった事がありまして、それほど苦労せずに思い出せました。
もしかしたら1024か2048辺りで割った方がいいのでしょうか？？（どっちにしろコヒーレンスは１以上にならないと思ったので気にしてませんでした）

フーリエスペクトルとパワースペクトルはなんとなくそんな感じがしてました。
ただ、なんとも自信が無くてルートを取ったりつけてたりして比較してました。

>アンチエリアジングフィルター
サンプリングエラーって奴ですよね。これは一応対策してます。

>コヒーレンスを求めるには平均化処理をしなければなりません。
「スペクトル解析」日野幹雄著　の本を買ってきたんですが、それでも「？？？」でした。（学が足りないのでほとんど理解できません）

まず、現状ですが、
信号AのＦＦＴ（Ar[i],Ai[i]) (実部,虚部)
信号BのＦＦＴ（Br[i],Bi[i])
パワースペクトルA[i] = Ar[i]^2 + Ai[i]^2
パワースペクトルB[i] = Br[i]^2 + Bi[i]^2
クロススペクトル[i] = (Ar[i] * Br[i] + Ai[i] * Bi[i])^2 + (Ar[i] * Bi[i] + Ai[i] * Br[i])^2
コヒーレンス[i]　＝　クロススペクトル[i]/(パワースペクトルＡ[i]　＊　パワースペクトルＢ[i]）
と言う式で、iを0～2047までループしてます。

多分これであってると思います。後は平均化処理だけじゃないかと思います。
で、平均化するのはコヒーレンス[i]だと遅すぎると思うのでパワースペクトルor(and)クロススペクトルor(and)フーリエスペクトルだと思います。（入力データを平均しても意味無さそうですし）

この平均なんですが、どれの平均ですか？？　そして、何個の平均を取れば良いのでしょうか？
例えば
Ar[i] = (Ar[i] + Ar[i+1] + ar[i+2]) / 3
のように、数個（３～１０個ぐらい？）の平均でやっていくのか、2048の平均をAr[0]～ar[2047]に入れてしまっていいのか、もっと良い平均のとり方があるのか、さっぱりです。

回答No.2

noname#65504

＃１です。

訂正補足です

「自己相関関数をフーリエ変換したもので」は自己相関関数ではなく相互相関関でした。

回答No.1

noname#65504

リンク先をじっくり読んだわけではないこと（リンク先に飛ばないものもあったので）、私自身は独学でかじった程度なので自信がないことを前提に呼んでください。

「スペクトル解析」日野幹雄著　朝倉書店 p63によると、コヒーレンスの定義は統一されておらず、２乗がつくものとそうでないものの２種類があり、文献を読む際にはどれを定義しているか確認する必要があるそうです。
だから、いろいろなサイトを読むよりも１つのサイトや書物に頼った方がよいと思います。

A)クロスパワースペクトル密度とクロススペクトル、クロスパワースペクトルは別物と言う事でしょうか？

　一般にクロスパワースペクトルと呼ばれるものは、クロスパワースペクトル密度あるいは、クロスパワースペクトル密度関数とも呼ばれています。
（「ディジタル信号処理の理論　１　基礎・システム・制御」谷萩隆嗣著　コロナ社　p113)
だからクロスパワースペクトルとクロススペクトルは別物、クロスパワースペクトルとは同じ物だと思います。

なお、パワーという用語は一般に２乗を示しています。

B)については記号の定義がよくわからないので省略します。

C)についてはディジタル信号上の処理の結果として、そういうことになっていると思います。

D)の計算は式の考えは正しいと思いますが、記号の付け方による大きな誤解があると思います。先ずクロススペクトルは自己相関関数をフーリエ変換したもので、フーリエ変換結果は虚数となりますので虚部を持ちます。
虚部と実部を同じ記号を使っているので、勘違いをしてしまったのだと思いますので、（a+bi)と(c+di)などと記号を変えて検討してみてください。自分の計算間違いに気づくと思います

先に述べたように同じ物を示すのに違う用語が使われたり、逆に同じ用語なのに違った定義がされていますので、いろいろな資料を当たるよりも１つの資料を参考にした方がよいと思います。部分部分をつまみ読みするよりも、ディジタル信号処理関係の本に書いてありますので１冊通して読んだ方がよいと思います。
ちなみに私はFTTアナライザーユーザーなので、ここをよく参照しています。
http://www.onosokki.co.jp/HP-WK/c_support/tech_term/cf_fft/cf3.htm

お礼

お礼が遅くなって申し訳ありませんでした。
>コヒーレンスの定義は統一されておらず
そうですか。参考になります。

>パワーという用語は一般に２乗を示しています。
知りませんでした。参考にさせていただきます。

>自分の計算間違いに気づくと思います
気づきました。相殺できませんでした。

>ディジタル信号処理関係の本に書いてありますので１冊通して読んだ方がよいと思います。
先日本屋に行ってきたのですが、コヒーレンスについて書かれている本が一冊しか見つかりませんでした。
音・振動のスペクトル解析 　工博 金井浩　ISBN：4-339-01105-3
本屋で読んだ感じではおそらくこの事について書いてあるんだと思いましたが、本に書いている事が理解できなかったためあきらめました。

で、ネットで調べた事とsemi-zzz様の解説を組み合わせ、あと少しのところまで来てるのではないか？と思っているのですが今悩んでます。
それは、ネット上のフリーソフトを使いコヒーレンスを求めたところ、同じサイン波同士だと全データが１になりました。おそらく正しいです。
自分でやる場合ですが、ＦＦＴを求めるとＳＩＮ波なのである周波数が1024で、他は０になります。
次に、パワースペクトルですが、ある周波数が1024で他は全部０でした。
クロススペクトルが自信ないのですが、計算式から行けばおそらくある周波数が1024＾2になり、他は全部０になると思います。
その後コヒーレンスを求めるのですが、分母（パワースペクトル）が０で計算できない上に分子も０なので０を出すしかないです。さらにある周波数は１どころか1024＾2
でも、0/0＝1と言うのもありなのかな？とも思ってます。（もしそれが可能なら完成に近づきますが）
で、小野測器のＨＰでは平均しなければ駄目と書いてましたし、他にも何かする事があるのかな？と思ってます。

とても参考になりました。
ありがとうございました。