最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

ウェルチ法を使用したストリーミング パワー スペクトル推定

この例では、ウェルチの平均修正ピリオドグラム法を使用して、時間領域入力のストリーミング パワー スペクトル推定を出力するブロックを紹介します。

はじめに

Spectrum Estimator ブロックには、選択肢としてフィルター バンク ベースのスペクトル推定とウェルチの平均修正ピリオドグラム法があります。

フィルター バンク法では、入力時間領域信号がフィルター バンクを使って異なる周波数ビンに分割され、各サブバンド信号の平均強度が計算されます。フィルター バンク ベースのスペクトル推定の詳細については、高分解能スペクトル解析を参照してください。

ウェルチ法では、入力の時間領域データがオーバーラップ可能なデータ セグメントに分割されます。各セグメントにウィンドウが適用され、次に、ウィンドウが適用されたシーケンスに基づき平均ピリオドグラムが計算されます [1]。データ セグメントの長さとウィンドウの選択により、推定の分解能帯域幅 (RBW) が決定されます。RBW は分解できる最小の正の周波数です。

この例で紹介したブロックは、ウェルチ アルゴリズムを使用していますが、このアルゴリズムは Spectrum Analyzer スコープ ブロックでも使用されています。Spectrum Analyzer はスペクトルの可視化のために使用します。この例のブロックは、推定されたスペクトルに直接アクセスする必要がある場合に (単に可視化するより) 役に立ちます。出力されたパワー スペクトルは、モデル内の他のブロックへの入力として使用されるか、あるいは事後処理のためにワークスペースに記録される場合があります。

ブロック ダイアログ

この例では、Spectrum Estimator ブロックを使用します。

このブロックには、以下のように 1 つの入力端子 (時間領域信号) と最大 5 つの出力端子があります。

1.推定パワー スペクトル P。

2.最大ホールド スペクトル P_max。各周波数ビンで、すべてのパワー スペクトル推定の最大値を維持することによって計算されます。この出力はオプションです。

3.最小ホールド スペクトル P_min。各周波数ビンで、すべてのパワー スペクトル推定の最小値を維持することによって計算されます。この出力はオプションです。

4.周波数のベクトル F。ここでパワー スペクトルが推定されます。この出力はオプションです。

5.ブロックの有効分解能帯域幅 RBW。この出力はオプションです。

ブロックのダイアログを以下に示します。このブロックのパラメーターは、Spectrum Analyzer ブロック上でアクセス可能なパラメーターと同じです。ブロック ダイアログにより、以下のようなスペクトル推定アルゴリズムのさまざまな機能を制御できます。

  • スペクトル タイプの選択 (パワーまたはパワー密度)

  • ブロックの分解能を制御する方法の指定 (RBW を直接指定するか、ウィンドウの長さを指定)

  • 推定器のウィンドウ関数の指定

  • 連続するセグメント間におけるオーバーラップ率の指定

  • 周波数範囲を、両側、片側または中心のいずれかとして選択

ブロック アーキテクチャ

このブロックは次の 2 つの段階で形成されています。

1.最初の段階では、入力は選択されたウィンドウの長さとオーバーラップ率に基づくセグメントにバッファーされます。この段階は、Buffer ブロックを使用して実装されます。

2.第 2 段階では、オーバーラップしたセグメントのパワー スペクトルが計算され、平均化されます。この段階は、dsp.SpectrumEstimator System object を使用して実装されます。

ウィンドウの長さの指定

以下に示す dspstreamingwelch モデルでは、44100 Hz でサンプリングされた、ノイズを含むチャープ信号のスペクトル推定に Welch Spectrum Estimator ブロックが使用されています。パワー スペクトル推定は、Array Plot スコープを使用して表示されます。スペクトルのピーク値と、ピークが発生する周波数が検出され、スコープ上に表示されます。推定の RBW も表示されます。さらに、比較および検証のために、Spectrum Analyzer スコープ ブロックも含まれています。

ブロックの周波数分解能の方法は、'Window Length (ウィンドウの長さ)' に設定されています。ウィンドウの長さは 1024 に設定されています。FFT 長、NFFT はウィンドウの長さと等価です。データは、サイドローブ減衰 60 dB のチェビシェフ ウィンドウを使用してウィンドウ処理されています。周波数範囲は片側です。この場合、スペクトル推定の長さは NFFT/2+1 = 513、計算区間は [0 Hz,22050 Hz] です。これに基づき、Array Plot スコープの "Sample increment" プロパティは、Fs/NFFT = 44100/(1024 * 1000) に設定されます。ここでは、周波数単位を KHz でスケーリングするため、インクリメントが 1000 で除算されます。該当する [コンフィギュレーション プロパティ] ウィンドウを開くことで、スコープの "Sample increment" プロパティにアクセスできます。

分解能帯域幅については、参考文献 [2] で以下のように与えられています。

ここで、N はウィンドウの長さ、enbw はウィンドウの等価ノイズ帯域幅、SL は選択したチェビシフ ウィンドウのサイドローブ減衰、Fs はサンプルレートです。この場合、RBW は 65.38 Hz になります。

モデルをシミュレートする際に、表示された RBW 値が Spectrum Analyzer スコープの下のバーに表示されている値と等価であることが確認できます。また、2 つのブロックのピーク測定値は同じになります。

非ゼロのオーバーラップの指定

前の節のモデルにはオーバーラップがありませんでした。dspstreamingwelch_overlap モデルでは、50% のオーバーラップがある Welch Estimation ブロックを使用します。モデルの他のパラメーターが前の節のパラメーターと同一であるため、RBW は変わらず 65.38 Hz です。

ウィンドウの長さが 1024、オーバーラップ率が 50% であることから、新しいデータ セグメントを形成するには 512 個の入力サンプルが必要になります。入力データの長さが 1024 であるため、新しいデータ フレームがそれぞれ 2 つの新しいピリオドグラムを生成し、ブロックの出力端子は入力端子の 2 倍の速さで動作します。

この場合、ウェルチ法の推定ブロックにはゼロ レイテンシがないことに注意してください。最初のスペクトル推定出力はバッファーの初期条件に基づき、eps と等価になります。したがって、Spectrum Analyzer スコープのスペクトルと測定値を一致させるため、Spectrum Analyzer の入力に Delay ブロックを挿入します。

Spectrum Analyzer の結果とウェルチ法の推定ブロックは、モデルをシミュレートすることで検証できます。

RBW の指定

dspstreamingwelch_rbw モデルでは、'Frequency Resolution Method (周波数分解能の方法)' パラメーターは RBW に設定されています。RBW Source は自動です。Spectrum Analyzer スコープ ブロックと同様に、このモードでは、分解能帯域幅は、指定した周波数範囲にわたって RBW 間隔が 1024 となるように選択されます。この例では範囲が 22050 Hz、RBW が 21.53 Hz になります。

データのバッファーに使用されるウィンドウの長さは、目的の RBW が得られるまで反復計算されます。この例でのウィンドウの長さは 3073 になります。この値を検証するには、このウィンドウの長さから結果として得られる RBW を以下のように計算します。

このモデルではハン ウィンドウが使用されていることに注意してください。

この場合、FFT 長 NFFT は奇数であり、3073 (ウィンドウ長) と等しくなります。周波数範囲が片側であるため、スペクトル推定は長さ (NFFT + 1)/2 で、計算区間は [0,44100/2) になります。Array Plot スコープの "Sample increment" プロパティは、Fs/NFFT = 44100/(3073 * 1000) KHz に設定されています。

ここでも、Spectrum Analyzer の結果とウェルチ法の推定ブロックは、モデルをシミュレートすることで検証できます。

参考文献

[1] 'Statistical Digital Signal Processing and Modeling', M.H. Hayes, Wiley, 1996.

[2] https://www.mathworks.com/help/dsp/ref/spectrumanalyzer.html