スペクトル密度推定

Signal processing technique

統計的信号処理において、スペクトル密度推定（SDE）または単にスペクトル推定と呼ばれる手法の目的は、信号の時間サンプル系列から信号のスペクトル密度（パワースペクトル密度とも呼ばれる）を推定することである。[ 1 ]直感的に^言えば、スペクトル密度は信号の周波数成分を特徴付ける。スペクトル密度を推定する目的の一つは、データ内の周期性を検出することであり、そのためには、これらの周期性に対応する周波数におけるピークを観測する必要がある。

一部のSDE手法では、信号が限られた数（通常は少数）の発生周波数とノイズで構成されていると仮定し、発生周波数の位置と強度を求めます。一方、成分の数について仮定を設けず、発生スペクトル全体を推定しようとする手法もあります。

概要

音声波形とその周波数スペクトルの例

周期的な波形（三角波）とその周波数スペクトル。220 Hzの「基本」周波数とそれに続く220 Hzの倍数（高調波）を示しています。

音楽のセグメントのパワースペクトル密度は、比較のために2つの異なる方法で推定されます。

スペクトル解析は、周波数領域解析またはスペクトル密度推定とも呼ばれ、複雑な信号をより単純な要素に分解する技術的なプロセスです。前述のように、多くの物理プロセスは、多数の個々の周波数成分の和として最もよく説明されます。周波数（または位相）に対する様々な量（振幅、電力、強度など）を定量化するプロセスはすべて、スペクトル解析と呼ばれます。

スペクトル解析は信号全体に対して行うことができます。あるいは、信号を短いセグメント（フレームと呼ばれることもあります）に分割し、個々のセグメントに対してスペクトル解析を行うこともできます。周期関数（など）は、この細分化に特に適しています。非周期関数を解析するための一般的な数学的手法は、フーリエ解析の範疇に入ります。 ${\displaystyle \sin(t)}$

関数のフーリエ変換は、元の信号に関するすべての情報を含む周波数スペクトルを生成しますが、形式は異なります。つまり、元の関数は逆フーリエ変換によって完全に再構成（合成）できるということです。完全な再構成のためには、スペクトルアナライザは各周波数成分の振幅と位相の両方を保存しなければなりません。これら2つの情報は、2次元ベクトル、複素数、または極座標における振幅と位相（つまり位相器）として表すことができます。信号処理でよく使われる手法は、振幅の2乗、つまりパワーを考えることです。この場合、結果として得られるプロットはパワースペクトルと呼ばれます。

フーリエ変換は可逆性を持つため、関数を時間ではなく周波数で表現したもの、つまり周波数領域表現と呼ばれます。時間領域で実行できる線形演算は、多くの場合、周波数領域でより簡単に実行できる同等の演算が存在します。周波数解析は、線形および非線形の両方を含む様々な時間領域演算の効果の理解と解釈を容易にします。例えば、周波数スペクトルに新しい周波数を作成できるのは、非線形または時間変動演算のみです。

実際には、周波数スペクトルを生成するほぼすべてのソフトウェアと電子デバイスは、信号のサンプルに対して動作し、完全な積分解の数学的近似値を提供する離散フーリエ変換(DFT) を使用します。 DFT は、ほぼ常に高速フーリエ変換(FFT) と呼ばれる効率的なアルゴリズムによって実装されます。 DFT の 2 乗振幅成分の配列は、ピリオドグラムと呼ばれるパワー スペクトルの一種であり、フィルタ インパルス応答やウィンドウ関数などのノイズのない関数の周波数特性を調べるために広く使用されます。 しかし、ピリオドグラムは、ノイズのような信号や、信号対雑音比の低い正弦波に適用された場合、処理ゲインを提供しません^{[なぜ? ] 。言い換えると、特定の周波数でのスペクトル推定値の分散は、計算に使用されるサンプル数が増えても減少しません。 これは、時間にわたって平均化 (}ウェルチ法^[2] ) するか、周波数にわたって平均化 (平滑化)することで軽減できます。ウェルチ法はスペクトル密度推定（SDE）に広く用いられています。しかし、ピリオドグラムに基づく手法では、一部の用途では許容できない小さなバイアスが生じます。そこで、次のセクションでは他の代替手法を紹介します。

テクニック

基本的なピリオドグラムの欠点を軽減するために、スペクトル推定のための他の多くの手法が開発されてきた。これらの手法は一般に、ノンパラメトリック法、 パラメトリック法、および最近ではセミパラメトリック法（スパース法とも呼ばれる）に分類できる。^[3]ノンパラメトリック手法は、プロセスが特定の構造を持っていると仮定せずに、プロセスの共分散またはスペクトルを明示的に推定する。基本的なアプリケーションで使用される最も一般的な推定値の一部（たとえば、ウェルチ法）は、ピリオドグラムに密接に関連したノンパラメトリック推定値である。対照的に、パラメトリック手法は、基礎となる定常確率過程が、少数のパラメータを使用して記述できる特定の構造（たとえば、自己回帰モデルまたは移動平均モデルを使用）を持っていると仮定する。これらの手法では、確率過程を記述するモデルのパラメータを推定することがタスクである。セミパラメトリック法を用いる場合、基礎となるプロセスはノンパラメトリックな枠組みを用いてモデル化され、モデルの非ゼロ成分の数が少ない（すなわち、モデルがスパースである）という仮定が追加されます。同様のアプローチは、欠損データの回復^[4]や信号再構成にも用いられます。

以下はスペクトル密度推定手法の一部です。

• 信号サンプルが時間的に不均等に間隔をあけられるノンパラメトリック法（記録が不完全な場合がある）
  • 既知の周波数に対する最小二乗法に基づく最小二乗スペクトル解析
  • ロンブ・スカーグル・ピリオドグラム、最小二乗スペクトル解析の近似
  • 非一様離散フーリエ変換
• 信号サンプルが時間的に均等に間隔をあけていなければならない（記録は完全でなければならない）ノンパラメトリック法：
  • ピリオドグラム、離散フーリエ変換の二乗係数
  • バートレット法は、信号の複数のセグメントから取得したピリオドグラムの平均をとり、スペクトル密度推定値の分散を減らす方法である。
  • ウェルチ法は、重複するセグメントを使用するバートレット法のウィンドウバージョンである。
  • マルチテーパーは、複数のテーパー、またはウィンドウを使用してスペクトル密度の独立した推定値を形成し、スペクトル密度推定値の分散を減らす、ピリオドグラムベースの方法です。
  • 特異スペクトル分析は、共分散行列の特異値分解を用いてスペクトル密度を推定するノンパラメトリック手法である。
  • 短時間フーリエ変換
  • 臨界フィルタは、ノイズ、不完全なデータ、機器応答関数を扱うことができる情報場理論に基づくノンパラメトリック手法である。
• パラメトリック手法（不完全なリスト）:
  • n番目のサンプルが前のp 個のサンプルと相関していると仮定する自己回帰モデル(AR) 推定。
  • 移動平均モデル(MA) 推定。n 番目のサンプルが前の p サンプルのノイズ項と相関していると想定します。
  • AR モデルと MA モデルを一般化する自己回帰移動平均(ARMA) 推定。
  • MUltiple SIGNAL CLASSIFIC (MUSIC) は、人気の高い超解像手法です。
  • 回転不変性技術による信号パラメータの推定(ESPRIT) は、もう 1 つの超解像法です。
  • 最大エントロピースペクトル推定は、鋭いピークなどの特異なスペクトル特性が予想される場合の SDE に役立つ全極法です。
• セミパラメトリック手法（不完全なリスト）
  • SParse反復共分散ベース推定（SPICE）推定^[3]およびより一般化された-SPICE ^{[5] 。} ${\displaystyle (r,q)}$
  • 反復適応アプローチ（IAA）推定。^[6]
  • Lassoは最小二乗スペクトル解析に似ていますが、スパース性を強制するペナルティが付きます。^[7]

パラメトリック推定

パラメトリックスペクトル推定では、信号が周波数とパラメータの関数であるスペクトル密度関数（SDF）を持つ定常プロセスによってモデル化されていると仮定します。^[8]すると、推定問題はこれらのパラメータを推定することになります。 ${\displaystyle S(f;a_{1},\ldots ,a_{p})}$ ${\displaystyle f}$ ${\displaystyle p}$ ${\displaystyle a_{1},\ldots ,a_{p}}$

パラメトリックSDF推定の最も一般的な形式は、次数の自己回帰モデル をモデルとして使用する。^{[8] : 392}ゼロ平均プロセスに従う信号シーケンスは、 次の式を満たす。 ${\displaystyle {\text{AR}}(p)}$ ${\displaystyle p}$ ${\displaystyle \{Y_{t}\}}$ ${\displaystyle {\text{AR}}(p)}$

$${\displaystyle Y_{t}=\phi _{1}Y_{t-1}+\phi _{2}Y_{t-2}+\cdots +\phi _{p}Y_{t-p}+\varepsilon _{t},}$$

ここで、は固定係数であり、は平均ゼロ、イノベーション分散を持つ白色ノイズ過程である。この過程のSDFは ${\displaystyle \phi _{1},\ldots ,\phi _{p}}$ ${\displaystyle \varepsilon _{t}}$ ${\displaystyle \sigma _{p}^{2}}$

$${\displaystyle S(f;\phi _{1},\ldots ,\phi _{p},\sigma _{p}^{2})={\frac {\sigma _{p}^{2}\Delta t}{\left|1-\sum _{k=1}^{p}\phi _{k}e^{-2\pi ifk\Delta t}\right|^{2}}},\qquad |f|<f_{N},}$$

サンプリング時間間隔とナイキスト周波数を指定します。 ${\displaystyle \Delta t}$ ${\displaystyle f_{N}}$

プロセスのパラメータとスペクトル密度を推定するアプローチはいくつかある： ^{[8] : 452-453} ${\displaystyle \phi _{1},\ldots ,\phi _{p},\sigma _{p}^{2}}$ ${\displaystyle {\text{AR}}(p)}$

• ユール・ウォーカー推定量は、ユール・ウォーカー方程式をプロセスに対して再帰的に解くことによって求められる。 ${\displaystyle {\text{AR}}(p)}$
• バーグ推定量は、ユール・ウォーカー方程式を通常の最小二乗問題の一種として扱うことによって求められる。バーグ推定量は、ユール・ウォーカー推定量よりも一般的に優れていると考えられている。^{[8] : 452}バーグは、これらを最大エントロピースペクトル推定と関連付けた。^[9]
• フォワード・バックワード最小二乗推定量は、プロセスを回帰問題として扱い、フォワード・バックワード法を用いてその問題を解きます。これはバーグ推定量と競合します。 ${\displaystyle {\text{AR}}(p)}$
• 最尤推定量は、最尤アプローチを用いてパラメータを推定します。これは非線形最適化を伴い、最初の3つよりも複雑です。

代替のパラメトリック手法には、移動平均モデル(MA) へのフィッティングや完全自己回帰移動平均モデル(ARMA) へのフィッティングが含まれます。

頻度推定

周波数推定とは、ノイズが存在する状況で、信号の成分の数に関する仮定に基づいて、信号の周波数、振幅、位相シフトを推定するプロセスです。 ^[10]これは、成分について事前に仮定を立てない上記の一般的な方法とは対照的です。

単音

最も大きな純音信号の周波数のみを推定したい場合は、ピッチ検出アルゴリズムを使用できます。

支配的な周波数が時間とともに変化する場合は、時間周波数表現で定義された瞬間周波数の推定が問題となる。瞬間周波数の推定法としては、ウィグナー・ヴィル分布や高次あいまい関数に基づく方法などがある。^[11]

受信信号（送信信号とノイズを含む）の（おそらく複雑な）周波数成分をすべて知りたい場合は、マルチトーンアプローチを使用します。

複数のトーン

信号の典型的なモデルは、白色ノイズが存在する場合の複素指数関数の和から構成され、 ${\displaystyle x(n)}$ ${\displaystyle p}$ ${\displaystyle w(n)}$

${\displaystyle x(n)=\sum _{k=1}^{p}A_{k}e^{in\omega _{k}}+w(n)}$。

のパワースペクトル密度は、ノイズによるスペクトル密度関数に加えて、インパルス関数で構成されます。 ${\displaystyle x(n)}$ ${\displaystyle p}$

周波数推定の最も一般的な方法は、ノイズ部分空間を特定してこれらの成分を抽出することです。これらの方法は、自己相関行列を信号部分空間とノイズ部分空間に固有分解することに基づいています。これらの部分空間が特定された後、周波数推定関数を用いてノイズ部分空間から成分周波数を求めます。ノイズ部分空間に基づく周波数推定の最も一般的な方法は、ピサレンコ法、多重信号分類法（MUSIC法）、固有ベクトル法、最小ノルム法です。

ピサレンコ法

${\displaystyle {\hat {P}}_{\text{PHD}}\left(e^{j\omega }\right)={\frac {1}{\left|\mathbf {e} ^{H}\mathbf {v} _{\text{min}}\right|^{2}}}}$

音楽

${\displaystyle {\hat {P}}_{\text{MU}}\left(e^{j\omega }\right)={\frac {1}{\sum _{i=p+1}^{M}\left|\mathbf {e} ^{H}\mathbf {v} _{i}\right|^{2}}}}$

固有ベクトル法

${\displaystyle {\hat {P}}_{\text{EV}}\left(e^{j\omega }\right)={\frac {1}{\sum _{i=p+1}^{M}{\frac {1}{\lambda _{i}}}\left|\mathbf {e} ^{H}\mathbf {v} _{i}\right|^{2}}}}$

最小ノルム法

${\displaystyle {\hat {P}}_{\text{MN}}\left(e^{j\omega }\right)={\frac {1}{\left|\mathbf {e} ^{H}\mathbf {a} \right|^{2}}};\ \mathbf {a} =\lambda \mathbf {P} _{n}\mathbf {u} _{1}}$

計算例

からまでの平均がゼロの時系列（離散時間）であると仮定します。これは有限個の周期成分（すべての周波数が正）の和であると仮定します。 ${\displaystyle x_{n}}$ ${\displaystyle n=0}$ ${\displaystyle N-1}$

$${\displaystyle {\begin{aligned}x_{n}&=\sum _{k}A_{k}\sin(2\pi \nu _{k}n+\phi _{k})\\&=\sum _{k}A_{k}\left[\sin(\phi _{k})\cos(2\pi \nu _{k}n)+\cos(\phi _{k})\sin(2\pi \nu _{k}n)\right]\\&=\sum _{k}\left[a_{k}\cos(2\pi \nu _{k}n)+b_{k}\sin(2\pi \nu _{k}n)\right]\end{aligned}}}$$どこ $${\displaystyle {\begin{aligned}a_{k}&=A_{k}\sin(\phi _{k}),&b_{k}&=A_{k}\cos(\phi _{k}).\end{aligned}}}$$

の分散は、上記のように平均ゼロの関数の場合、次のように表される。 ${\displaystyle x_{n}}$

$${\displaystyle {\frac {1}{N}}\sum _{n=0}^{N-1}x_{n}^{2}.}$$

これらのデータが電気信号から取得されたサンプルである場合、これは平均電力になります (電力は単位時間あたりのエネルギーであるため、エネルギーが振幅の 2 乗に類似しているとすれば、これは分散に類似しています)。

ここで、簡単にするために、信号が時間的に無限に拡張すると仮定すると、次の限界に達します。平均電力が制限されている場合 (実際にはほとんどの場合に当てはまります)、次の限界が存在し、データの分散となります。 ${\displaystyle N\to \infty .}$

$${\displaystyle \lim _{N\to \infty }{\frac {1}{N}}\sum _{n=0}^{N-1}x_{n}^{2}.}$$

ここでも、話を簡単にするために連続時間を想定し、信号が両方向に無限に広がると仮定する。すると、これらの2つの式は

$${\displaystyle x(t)=\sum _{k}A_{k}\sin \left(2\pi \nu _{k}t+\phi _{k}\right)}$$

そして

$${\displaystyle \lim _{T\to \infty }{\frac {1}{2T}}\int _{-T}^{T}x(t)^{2}\,dt.}$$

の二乗平均平方根はなので、 の分散はである。したがって、の周波数成分からの平均電力への寄与はである。これらすべての寄与を合計すると の平均電力となる。 ${\displaystyle \sin }$ ${\displaystyle 1/{\sqrt {2}}}$ ${\displaystyle A_{k}\sin(2\pi \nu _{k}t+\phi _{k})}$ ${\displaystyle {\tfrac {1}{2}}A_{k}^{2}.}$ ${\displaystyle x(t)}$ ${\displaystyle \nu _{k}}$ ${\displaystyle {\tfrac {1}{2}}A_{k}^{2}.}$ ${\displaystyle x(t).}$

すると、頻度の関数としての電力は、その統計的累積分布関数は、 ${\displaystyle {\tfrac {1}{2}}A_{k}^{2},}$ ${\displaystyle S(\nu )}$

$${\displaystyle S(\nu )={\frac {1}{2}}\sum _{k:\nu _{k}<\nu }A_{k}^{2}.}$$

${\displaystyle S}$はステップ関数であり、単調に減少しない。そのジャンプは の周期成分の周波数で発生し、各ジャンプの値は当該成分のべき乗または分散となる。 ${\displaystyle x}$

分散は、データ自身との共分散です。ここで、同じデータについて、ラグ を加えた場合、との共分散を求め、これを信号（またはデータ） の自己相関関数として定義することができます。 ${\displaystyle \tau }$ ${\displaystyle x(t)}$ ${\displaystyle x(t+\tau )}$ ${\displaystyle c}$ ${\displaystyle x}$

$${\displaystyle c(\tau )=\lim _{T\to \infty }{\frac {1}{2T}}\int _{-T}^{T}x(t)\,x(t+\tau )\,dt.}$$

それが存在する場合、それは の偶関数です。平均検出力が有界である場合、 はどこにでも存在し、 は有限であり、 によって有界となります。これは、データの平均検出力または分散です。 ${\displaystyle \tau .}$ ${\displaystyle c}$ ${\displaystyle c(0)}$

は と同じ周期を持つ周期成分に分解できることが示されます。 ${\displaystyle c}$ ${\displaystyle x}$

$${\displaystyle c(\tau )={\tfrac {1}{2}}\sum _{k}A_{k}^{2}\cos(2\pi \nu _{k}\tau ).}$$

これは実際には、異なる周波数にわたる のスペクトル分解であり、 の周波数にわたる電力の分布に関連しています。 の周波数成分の振幅は、信号の平均電力への寄与です。 ${\displaystyle c}$ ${\displaystyle x}$ ${\displaystyle c}$

この例のパワースペクトルは連続していないため、微分係数を持たず、したがってこの信号はパワースペクトル密度関数を持ちません。一般に、パワースペクトルは2つの部分の和になります。1つは、この例のように連続しておらず密度関数を持たない線スペクトル、もう1つは完全に連続しており密度関数を持つ残差スペクトルです。

参照

• 多次元スペクトル推定
• ピリオドグラム
• シグスペック
• スペクトログラム
• 時間周波数解析
• 時間周波数表現
• ウィットル尤度
• スペクトルパワー分布

参考文献

1. P Stoicaと R Moses、「信号のスペクトル分析」、Prentice Hall、2005 年。
2. Welch, PD (1967)、「パワースペクトルの推定における高速フーリエ変換の利用：短い修正ピリオドグラムの時間平均に基づく手法」、IEEE Transactions on Audio and Electroacoustics、AU-15 (2): 70– 73、Bibcode :1967ITAE...15...70W、doi :10.1109/TAU.1967.1161901、S2CID  13900622
3. Stoica, Petre; Babu, Prabhu; Li, Jian (2011年1月). 「分離可能モデルにおけるスパースパラメータ推定の新手法と不規則サンプリングデータのスペクトル解析への応用」IEEE Transactions on Signal Processing . 59 (1): 35– 47. Bibcode :2011ITSP...59...35S. doi :10.1109/TSP.2010.2086452. ISSN  1053-587X. S2CID  15936187.
4. Stoica, Petre; Li, Jian; Ling, Jun; Cheng, Yubo (2009年4月). 「非パラメトリック反復適応アプローチによる欠損データ回復」. 2009 IEEE International Conference on Acoustics, Speech and Signal Processing . IEEE. pp.  3369– 3372. doi :10.1109/icassp.2009.4960347. ISBN 978-1-4244-2353-8。
5. Sward, Johan; Adalbjornsson, Stefan Ingi; Jakobsson, Andreas (2017年3月). 「スパース反復共分散ベース推定量の一般化」. 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) . IEEE. pp.  3954– 3958. doi :10.1109/icassp.2017.7952898. ISBN 978-1-5090-4117-6. S2CID  5640068。
6. Yardibi, Tarik; Li, Jian; Stoica, Petre; Xue, Ming; Baggeroer, Arthur B. (2010年1月). 「音源定位とセンシング：加重最小二乗法に基づくノンパラメトリック反復適応アプローチ」. IEEE Transactions on Aerospace and Electronic Systems . 46 (1): 425– 443. Bibcode :2010ITAES..46..425Y. doi :10.1109/TAES.2010.5417172. hdl : 1721.1/59588 . ISSN  0018-9251. S2CID  18834345.
7. Panahi, Ashkan; Viberg, Mats (2011年2月). 「LASSOベースのDOA推定法の解像度について」. 2011 International ITG Workshop on Smart Antennas. IEEE. pp.  1– 5. doi :10.1109/wsa.2011.5741938. ISBN 978-1-61284-075-8. S2CID  7013162。
8. パーシバル, ドナルド B.; ウォルデン, アンドリュー T. (1992). 『物理応用のためのスペクトル解析』ケンブリッジ大学出版局. ISBN 9780521435413。
9. Burg, JP (1967)「最大エントロピースペクトル解析」、探査地球物理学会第37回会議議事録、オクラホマシティ、オクラホマ州。
10. Hayes, Monson H.,統計的デジタル信号処理とモデリング, John Wiley & Sons, Inc., 1996. ISBN 0-471-59431-8。
11. Lerga, Jonatan. 「信号瞬時周波数推定法の概要」（PDF） . リエカ大学. 2014年3月22日閲覧。

さらに読む

• ポラット、B. (1994). 『ランダム信号のデジタル処理：理論と方法』 プレンティス・ホール. ISBN 978-0-13-063751-2。
• プリーストリー、MB（1991）『スペクトル解析と時系列』アカデミック・プレス、ISBN 978-0-12-564922-3。
• ストイカ, P.; モーゼス, R. (2005).信号のスペクトル解析. プレンティス・ホール. ISBN 978-0-13-113956-5。
• Thomson, DJ (1982). 「スペクトル推定と高調波解析」. Proceedings of the IEEE . 70 (9): 1055–1096 . Bibcode :1982IEEEP..70.1055T. CiteSeerX  10.1.1.471.1278 . doi :10.1109/PROC.1982.12433. S2CID  290772.

Retrieved from "https://en.wikipedia.org/w/index.php?title=Spectral_density_estimation&oldid=1303898846"