![]() |
|||||||||||||||
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
||||||
当事業所では海洋・水産・気象関連のデータ処理・解析や数値計算を主に行っていますが、他分野の諸計算・統計解析業務についてもお受け致します。個人やNPO、NGOのお客様も歓迎致します。
お問合せは |
![]() |
||||||||||||||
パワースペクトル密度計算 (スペクトル解析) 図を含む本文書の著作権はCygnus Research Internationalに帰属します。無断転載はお断り致します。 PD001A/B ユーザーガイド 1/7 目次 パワースペクトル密度関数(以下PSDと略)は、データ中に周期的な変動が存在する場合にその周期(あるいは周波数)や振幅を調べるためによく用いられます。例として図1-1(a)及び(b)に海洋の流速データの時系列及びそれより計算したPSDを示します。PSDの図(図1-1(b))では1日周期(図中でDと表示)及び半日周期(同SD)の潮汐成分の存在が容易にわかりますが、時系列の図(図1-1(a))ではこれらの変動があるかどうかはまったくわかりません。図1-1(c)は図1-1(a)の最初の240データのみを取り出して拡大した図で、この期間中に1日周期の潮汐は24回、半日周期の潮汐は48回あることになるのですが、この図でも両潮汐成分は明確ではありません。 本小額定型計算業務PD001シリーズではお客様のデータの数値の変動を異なった周波数の三角関数に分解し、それらの振幅を計算することによってPSDを求めるという一般によく使われる方法を採用しています。図1-2に三角関数を示しますが、もしお客様のデータの変動に鋭い角、刺(スパイク)や階段状(矩形波を含む)といった変動がある場合は、それらの変動がたとえ非常に正確な周期を持っていても本業務で計算した結果の解釈にはやや注意が必要となる可能性があります。->ここ<-で色々なデータのPSD例を記述していますので、よろしければ御参照下さい。 PSDを計算する際には事前に決めなければならない項目が幾つかありますが、大抵の場合は実際に計算を行ってその結果を調べてみなければそれらの設定が適切かどうかよくわからないのが現実です。そこでPD001Aでは下の表1Aに示すように1注文で9ケースの異なった設定で計算を行い、お客様の側でそれらの結果を比べた上でその中から最適な結果を選択して頂けるようにしています。ただし、PD001Bでは表1Bに示すように1注文で1ケーズの計算のみを行います。 以下では1章では価格を含む本業務の要約を、2章では計算結果を格納したファイルや図といった本業務のプロダクツの内容の詳細及びその利用例を、3章ではマイクロソフト社のエクセルを用いた入力データファイルの作製方法を、4章では本業務でお客様にご指定頂く各種設定についての解説を記述しています。1章の内容は簡潔なものとしていますが、詳細は2章以後に詳しく記述していますので、1章の内容でご不明な点があった場合は2章以後もお読み頂きますようお願い致します。 本文書では周波数と周期を適宜入れ換えて使う場合がありますが、周期は周波数の逆数(周期=1/周波数)で、例えば周波数が2サイクル/秒の場合の周期は1/2=0.5秒となり、周波数が高いと周期は短くなり、逆に周波数が低いと周期は長くなります。また、本文書では時間領域と周波数領域という言葉を用いていますが、これは図1-1のような場合に水平軸が時間なのか周波数なのかといったような違いで、図1-1(a)は時間領域でのデータの表示で、図1-1(b)はそれに対応する周波数領域でのデータの表示となります。本文書では時系列データを扱っていますが、もしお客様のデータが、例えば物質の色を直線に沿ってスキャンしたような何かの1次元の空間分布である場合は時間を空間に置き換えて下さいますようお願いします。 本文書はPSDの解説を目的としたものではなく、当時業所が提供する小額定型計算業務PD001シリーズの解説を目的としたものですので、数式や専門用語の使用はできるだけ避け、PSDの数学的に詳細な定義や導き方といった事よりは実際の計算結果がどのようなものになるのか、また実際の計算精度はどの程度かといった記述が中心となっています。また、一部の表現は数学や統計学的には厳密ではないことをあらかじめお断り致します。 (1章) PD001シリーズの要約 プロダクツはお客様にお渡しする計算結果です。なお、PD001Aでは位相についてはケースA1でのみ、振幅についてはケースA1、B1及びC1でのみ計算致します。これは位相や振幅については他のケースでは計算する意味が無いためで、その詳細については後述致します。PD001Bをご発注の場合も位相や振幅に関してはこれに準じます。計算結果を格納したファイルの各項目については2章に記述していますので、そちらをご参照お願い致します。なお、図についてはPD001Aをご発注の場合は各ケースのPSDの図をPDF形式でご提供致しますが、PD001Bをご発注の場合は図はご提供致しません。表1CにPD001AとPD001Bの差異を示します。 PD001シリーズで計算するPSDは片側PSDとなりますが、 片側PSDについては(2-1-2-3)に記述しています。 PSDが計算される周波数は、最低が0、最高が1/(2xサンプリング間隔)となり、各周波数間の間隔は一定で1/(サンプリング間隔xデータ数)となり、データ数/2+1点の周波数となります。なお、ここでサンプリング間隔というのはデータが取得された時間間隔です。例として、もしサンプリング間隔が5秒、データ数が200の場合は、PSDが計算される周波数は0サイクル/秒(0Hz)から1/(2x5)=1/10=0.1サイクル/秒(0.1Hz)でその間の周波数間隔は1/(5x200)=1/1000=0.001サイクル/秒(0.001Hz)となり、総計200/1+1=101周波数となります。 (1-2) お客様がご指定できる計算設定項目の要約 (1-2-2) PSDの値の信頼区間のパーセンテージ; 全ケースで同一の値を指定 (1-2-3) テューキー窓関数のテーパー比; ケースB1、B2及びB3(PD001A)のみ適用 値は個別に指定 (1-2-4) 周波数領域での平滑化(以下FDSと略)ビン巾; ケースA1、B1及びC1の3ケースを除く全ケースで個別に指定 (1-3) 計算手法 (1-4) 入力データ(お客様のデータ)について その他の条件としては、データのサンプリング間隔は一定でなければなりません。サンプリング間隔が一定ではない場合は当事業所で適当な内挿方法を用いてサンプリング間隔を一定としたデータを作成することも可能ですが、その場合は別途料金が加算となります。なお、データの個数については100万程度でも処理可能ですが、奇数個の場合はこちらで値が0のデータを1個付加し偶数個にしてから処理を行います。この方法はゼロパディングと呼ばれ、計算の精度等には影響を与えません。この作業については追加料金は発生致しません。 PD001シリーズではお客様のデータのチェックは基本的には行いません。したがいまして、御発注前にお客様の側でデータが適切なものかどうかを図1-1(a)のような時系列図を作成してチェックされることをお勧めいたします。もしお客様のデータに異常な値が含まれている場合、当事業所でそれらの異常な値を検出しそれらの値を無効データとして周辺の値より内挿して置き換えるような処理も可能ですが、その場合は異常な値の明確な定義が必要となり、また別途加算料金が発生いたします。なお、データの数値の単位をご連絡頂いた場合はPD001Aで作成する図のラベル等に使用致します。単位が不明な場合は図のラベルの単位は無記入となります。 御見積もりは無料ですので、データファイルがここや3章でご説明するような形式でない場合や異常な値の検出その他の前処理が必要な場合の追加料金については事前にお問い合わせ頂きますようお願い致します。 (1-5) プロダクツ (1-5-1) 結果ファイル (1-5-2) PDF形式の図 (1-6) 料金、納期およびご発注手順 PD001Bの料金は1データセットにつき1200円(税込み、Paypal社手数料込み)となります。もし複数のデータセットの処理をご要望の場合は2セット目より1セット当たり600円(税込み)加算となります。なお、複数データセットの処理料金に対する割引は同時にご発注頂いた場合のみ適用いたします。 納期は処理するデータセットの数(個々のデータセットに含まれるデータ数ではありません)にもよりますが、通常はご発注確定後1-2営業日程度となります。なお、ご発注確定は以下の手順では(手順4)に該当致します。 ご発注手順は以下のとおりとなります。 データの送受をメール添付以外の方法で行うことは可能ですが、別途作業料金および送料等の実費が加算となります。 (1-7) 個人情報、お客様のデータに関する情報保護 お客様のデータおよびPD001シリーズの計算結果は業務完了後1週間程度で削除致しますので、PD001シリーズの計算結果はお客様の側で保管頂きますようお願い致します。なお、業務完了日時はお客様のお支払いが確認できた日時となります。業務管理のためお客様のメールアドレス、ご発注日時、業務完了日時及びPD001シリーズの制御用ファイルは保管致しますが、この制御用ファイルにはお客様のデータは一切含まれません。 (2章) プロダクツの詳細及びその使用例 (2-1) 結果ファイル(本業務で納入するファイル) (2-1-1) 結果ファイルのフォーマット(形式) (2-1-2) 結果ファイルに記述する計算結果の内容 (2-1-2-1) 周波数 (Frequency) 周波数の単位はサンプリング間隔の単位によって決まります。例えば、サンプリング間隔が2日である場合の周波数の単位はサイクル/日です。周波数の単位をサンプリング間隔の単位と異なるものにしたい場合(前の例で周波数の単位をサイクル/秒に変更等)は発注時に明記下さいますようお願い致します。なお、サンプリング間隔をご連絡頂かない場合はサンプリング間隔は単位無しの1.0とみなし、周波数の単位はサイクル/単位無しの1.0となります。 一般にPD001シリーズのように一定周波数間隔でPSDを計算する場合、計算されたPSDはPSDが計算された周波数を中心に一定の周波数巾を持つ領域のPSDを表わすものという考え方をする事が多くあります。この周波数領域は通常ビンと呼ばれ、ビンの周波数巾はPSDの周波数間隔と同じです。上の例を用いると3番目のビンの中心周波数は0+((3-1)x周波数間隔)=2x0.001=0.002サイクル/秒で、このビンの周波数の下限は中心周波数 ー(周波数間隔/2)=0.002-(0.001/2)=0.002-0.0005=0.0015サイクル/秒で、上限は中心周波数 +(周波数間隔/2)=0.002+0.0005=0.0025サイクル/秒となります。本ドキュメントではこのビンの概念を所々で使用致します。ただし、ここで1つ注意しなければならない点はPSDや振幅スペクトルそのものはあくまでもこの結果ファイルに書き出される周波数において計算された値であって、各ビンに含まれる周波数域の平均のようなものではありません。 (2-1-2-2) 周期 (Period) (2-1-2-3) PSD (PSD) (2-1-2-4) PSDの信頼区間 (Conf Int) (2-1-2-5) 分散のパーセンテージ (% of Var) あるビンのPSDにビンの周波数巾(ビンによらず一定の値です)を掛けた値はそのビンに含まれる分散になり、周波数ゼロを除くすべてのビンの分散を足し合わせた値は平均値を除去したデータの通常の分散と同じになります。これを数式で表現すると この式により、あるビンのPSDにビンの周波数巾をかけた値を時系列データの分散で割った値は、そのビン内の変動が全体の変動のどの程度を含むかを示します。これにさらに100を掛けたものが分散のパーセンテージで、これを数式で表現すると (2-1-2-6) 振幅 (Amplitude) PSD計算時に窓関数を適用した場合は計算された振幅に対し補正を行いますが、(4-2-8)の記述のように振幅に対する補正はPSDに対して行う補正と少し異なります。このため窓関数を適用しない場合については振幅はPSDの2倍に周波数巾をかけたものの平方根と等しくなりますが、窓関数を適用した場合はその関係が成り立ちません。なお、振幅に対してFDSを適用すると無意味な結果にしかなりませんので、表1Aや表1Bに示すようにFDSを適用したケースについては振幅は結果ファイルには書き出しませんが、これについての詳細は(4-3-2)に記述しています。 (2-1-2-7) 位相 (Phase) (2-1-3) 結果ファイルに書き出された情報の利用方法及び使用例 ここでは、すべてデフォルト設定でPD001Aによって計算を行い、その結果よりFDSは無適用でハニング窓関数を適用した場合のPSDの図を作成したいと仮定します。最初に行わなければならない作業は作図を行うソフトウェアにPD001Aの結果ファイルを読込ませる作業です。エクセルをご使用になる場合は、WindowsのエクスプローラやMacのファインダの画面で目的とする結果ファイルをダブルクリックして頂ければ、通常はエクセルが起動して自動的に読込まれます。ただし、ファイルの拡張子"csv"をエクセル以外に関連づけされている場合はまずエクセルを起動してから目的の結果ファイルを開いて頂くことになります。エクセル以外のソフトウェアをご使用になる場合はあらかじめお好みのテキストエディタで結果ファイルを読込み、最初の行(ヘッダ行)を削除したファイルを別途作成したほうが便利かもしれません。ただし、当事業所ではお客様の結果ファイルの保存は行いませんので、結果ファイルの内容を変更される場合は事前にコピーを作成してお客様ご自身で保存されることをお勧めいたします。 次に作図する範囲を周波数または周期で決め、その範囲がどの行からどの行になるかを結果ファイルの第1列(周波数)または第2列(周期)を参照して探します。エクセルではA列またはB列になります。次にこの例での作図対象はFDSは無適用でハニング窓関数を適用した場合のPSDですが、それがどの列にあたるのかを表2Aを用いて探します。この作図対象はケースC1のPSDに相当し、表2Aより30列目(エクセルではAD)になることがわかります。後はこの30列目で作図したい範囲の行を選択し、エクセルその他のソフトウェアで作図するだけです。 なお、作図上の注意事項ですが、(2-2)に記述するようにPSDの図を作成する場合は一般的には縦軸に対数スケールを用います。これは通常はPSDの変化の幅が非常に大きく、一般的によく使う均等(線形)スケールで作図すると非常に大きな値以外は見えなくなってしまうためです。また、図1-1(b)に示すように比較的広い周波数幅のPSDの作図を行うと高周波域でのPSDの大きな変動が目立ち過ぎる場合があります。このような場合は低い周波数域ではケースA1、B1またはC1の結果を用い、中間の周波数域ではケースA2、B2またはC2の結果を用い、高い周波数域ではケースA3、B3またはC3の結果を用いて作図すると見やすい図ができる場合があります。図2-1はこのような方法を用いて作成した図で、図中の黒の縦線はケースを切り替えた境界を示します。図2-1中の信頼区間については次の(2-1-3-2)に記述致します。この図のPSDの計算に用いたデータは図1-1(b)の計算に用いたデータと同じですが、図1-1(b)では全周波数域においてFDSビン巾7のFDSを適用した結果となっています。なお、図2-1のような図を作成したい場合はあらかじめ結果ファイルよりそれぞれのケースから作図に使用したい範囲の値だけをコピーして別のワークシート等にペーストして作図したい値だけを含むような新しいワークシート等を作成したほうが便利です。 (2-1-3-2) PSD及びPSDの信頼区間の使用方法 特定の変動の周波数を出来るだけ正確に知りたい場合はFDSを適用しないため周波数解像度が高いケースA1またはB1のどちらかを使用されればよいかと思われます。一般的にはケースA1やB1の結果には非常に沢山のピーク(変動が大きい周波数)が現れますので、目的のピークを確定するためにはケースA2やB2の結果も参照されたほうがよいかもしれません。ケースB1の解像度はテューキー窓関数のテーパー比で変わりますが、デフォルト設定の場合はケースA1の解像度とほとんど同じで、ケースC1は(4-2-6)で記述するようにハニング窓関数が持つFDSに似た効果により周波数解像度がケースA1及びB1に比べやや低下します。なお、FDSを適用した場合のPSD(デフォルト設定の場合はケースA1、B1及びC1を除く全ケース)は(4-3-1)で記述するように値が小さくなりますのでご注意が必要です。 ここでケースA1についてですが、統計の専門家は一般的にはケースA1の結果は使用しません。その理由を簡単にご説明いたしますと、まず、統計の専門家は通常はお客様が実際にお持ちのデータは統計的性格が同じような無数のデータセットの中の単なる1つのデータセットと見なし、お客様のデータを用いてPD001シリーズで計算したPSDはその無数のデータセットによってあらわされる真のPSDの単なる推定の1つにしか過ぎないと考えます。そうすると、計算されたPSDは推定された値にしかすぎないのですから、真のPSDに何らかの誤差が加わった値であると考えられます。なお、この誤差というのは補足1に記述するような数値計算を行うと常に生じる計算上の誤差とは無関係です。 ここで、真のPSDは通常は手持ちのデータから計算できるものではなく、場合によれば仮想的なものです。統計の専門家はできるだけ精確に真のPSDを知ることを望むわけですが、統計の理論を導入するとケースA1の結果に含まれていると推定される誤差は非常に大きく、従ってケースA1の結果は使用できないという結論に至ります。ケースB1の結果についてもテーパー比がデフォルトの値であれば、統計の専門家は多分使用しないでしょう。他のケースの結果については、統計の理論によれば誤差を減少させる効果がFDSにあるので使用しても差し支えがないと考えられます。ここでこの誤差の大きさですが、これは計算されたPSDの確かさの範囲という形でPSDの信頼区間によって示されます。なお、こういった議論や実際のPSDの信頼区間の計算は基本的にはデータの単純な平均によって得られるサンプル平均値(推定された平均値)とデータが属する母集団とよばれるものの平均(真の平均値)との関係と似ていて->ここ<-にもうすこし詳しく記述しています。(実際は母集団の平均ではなく分散の推定と同じです) ここで、もしお客様が実際には計算できない言わば仮想的な真のPSDにはまるで興味がなく、お持ちのデータのPSDが知りたいという場合は上のような統計の専門家の議論は無意味と言えるかもしれません。したがってそのような場合はケースA1の結果を使用してもかまわないと思われますが、その際にはPSDの信頼区間は意味をもちません。これはPD001シリーズで計算するPSDの信頼区間は上に記述したような統計の専門家の考えを基礎にしているからです。また、(2-1-3-3)で使用しているようなデータについては、そもそもデータの値に統計的にしか記述できない要素がまったく無く、また、データ長やデータの開始時間も決められているため上に記述したような誤差がデータには一切含まれていません。したがってこのような場合についてもケースA1の結果を使うことにはまるで問題がありません。なお、このような場合はPD001シリーズで計算するようなPSDの信頼区間は無意味なものとなります。これに対し(2-1-3-5)で使用しているようなデータの場合は、本来は統計的にしか記述できない要素である式(7)中のノイズが含まれているため、ケースA1の結果を使用される場合は上に記述したような点について留意されたほうがよいと思われます。また、(2-1-3-3)で使用しているようなデータについてもデータの開始時間やデータ長が不規則に変わるような場合は統計的にしか記述できない要素が加わります。 PSDの信頼区間の使用方法については具体例を用いて説明致します。図2-2(a)は細かいピークが非常に多いが全体としては比較的平坦なPSDの例を示します。この図中のAは最大のピークでBは2番目に大きなピークです。ここで問題となるのはこれらのピークがはたして特別なのか、それとも他の無数のピークと同様のものなのかという点ですが、それを判断するためには通常はAやBのピークが他と比べて特別に高いかどうかを調べます。そこで、最初に行う事は比較の基準となる”普通の高さ”を決めることになりますが、ここで、”普通の高さ”の値を”底”と呼ぶ事にします。この”底”の決め方はいろいろあるのですが、ここでは全周波数域でのPSDの95%の信頼区間の上限の平均値を”底”とします。なお、他の”底”の決め方としてはPSDそのものの平均、あるいはPSDそのものの平均にPSDの標準偏差の3倍を加えたもの等が考えられます。またPSDが図2-2(a)のように平坦ではなく傾斜が明確に見える場合等は大雑把に傾斜を求めてそれを”底”にする、あるいは調べたいピーク周辺にのみ平均を計算する周波数域を限定する等も考えられます。その他、周波数による特定の分布を仮定し、実際のPSDの計算結果を用いて調整したものを”底”にするといったことも行われます。いずれの方法を採用するにしろ、”底”を決める場合はその方法が妥当であることが大事ですが、そのためには最初に図2-2(a)のような図を一度眺めてみる必要があります。 さて、ここで図2-2(b)にPSD(青線)、PSDの95%の信頼区間の上限(赤線)、PSDの95%の信頼区間の下限(紫線)及び”底”(水平の赤線)をA周辺のみ拡大して示します。この図によれば、Aの信頼区間の下限(紫線)は”底”(水平の赤線)より上にありますが、このような場合は95%の信頼限界でAは”普通”に比べて大きい、すなわち特別と見なせます。これに対し図2-2(c)はB周辺を拡大表示したものですが、こちらの場合はPSD(青線)及びその信頼区間の上限(赤線)は”底”より高いのですが、信頼区間の下限は”底”より低くなっています。このような場合は95%の信頼限界でBは特別とは見なせないということになります。図2-2(a)によれば他の無数のピークはすべてBより低いのですから、このデータではAのみ特別で他はすべて普通となります。 上の例ではPSDは比較的平坦でしたので縦軸に均等スケールを採用しましたが、実際にはPSDの変動巾が大きいため縦軸を対数スケールにしたほうがよい場合がしばしばあります。縦軸を対数スケールにすると図中での見た目のPSDの信頼区間の上限と下限の巾はFDSの巾を一定に保つ限りどこでも同じになります。したがって、図1-1(b)の場合では信頼区間は図中の一カ所にだけ表示すればよく、図2-2のように全周波数域に折れ線で表示する必要はありません。図2-1の場合は周波数域によって異なるFDSの結果を使用しているため、計3カ所に信頼区間の巾を表示しています。 さて、ここで図2-3にこのような場合の信頼区間の使用例を示します。この例ではPSDは対数スケールを採用しても変動が大きく、全周波数域で”底”を決めるには曲線を採用しなければいけませんが、ここでは短い赤の水平の線で囲まれたピークを調べるためにかなり大雑把にその麓にあたる周辺の値を”底”としています。この図の例ではこの赤の水平線で囲まれたピークの高さは信頼区間の巾より高いので、特別と見なせますが、他のピークはそれぞれの麓を”底”にするとどれも高さが信頼区間の巾より低いので特別とは見なせません。 なお、信頼区間による判断は色々な仮定に基づいた一般的な統計理論の適用にすぎず、データの変動を実際に生じさせている物理的その他の過程は一切考慮されていません。従って、信頼区間を用いて特別と判断されたピークでもそれが実際に意味のあるものという保証はできませんし、また、信頼区間を用いて普通と判断されたピークが実際に意味のないノイズのようなものという断定もできません。以上により、信頼区間のみでPSDの図上にあらわれたピークの重要性を判断するのはあまりお勧めできません。PSDの信頼区間はPSDの図上に現れた無数のピークの中から特定のピークを選んだ場合において、その特定のピークを選んだ事を正当化するために使用する程度に止めておいたほうが無難ではないかと思われます。 |
|||||||||||||||