MATLAB でのパワー スペクトルの推定
時間領域信号のパワー スペクトル (PS) とは、有限長のデータをベースに、信号内に含まれているパワーを周波数上に分布したものです。信号の周波数領域表現は、時間領域表現よりも解析が容易である場合が多くあります。ノイズ キャンセリングやシステム同定といった多くの信号処理アプリケーションは、周波数に特化した信号修正に基づいています。パワー スペクトル推定の目的は、時間サンプルのシーケンスから信号のパワー スペクトルを推定することです。信号に関する既知の情報に応じて、パラメトリック アプローチまたはノンパラメトリック アプローチが推定法として使用され、時間領域解析または周波数領域解析に基づいて推定が行われます。たとえば、一般的なパラメトリック手法では、観測値を自己回帰モデルに適合させます。一般的なノンパラメトリック手法はピリオドグラムです。パワー スペクトルは、ウェルチ法やフィルター バンク法といったフーリエ変換法を使用して推定されます。比較的長さが短い信号の場合、フィルター バンク アプローチによって、より高い分解能とより正確なノイズ フロアをもつスペクトル推定が生成され、ウェルチ法よりも正確なピークを示します。スペクトル漏れは低いかまったくありません。こうした利点には、計算量の増加と追跡の低速化という代償が伴います。これらの方法の詳細については、Spectral Analysisを参照してください。最大エントロピー法のようなその他の手法を使用することもできます。
MATLAB® では、spectrumAnalyzer
オブジェクトを使用して動的な信号のリアルタイム スペクトル解析を実行することができます。オブジェクト関数 isNewDataReady
および getSpectrumData
を使用して、スペクトル データをスペクトル アナライザーに表示し、そのデータをワークスペース変数に格納することができます。あるいは、dsp.SpectrumEstimator
System object™ を使用した後で dsp.ArrayPlot
オブジェクトを使用することでスペクトル データを表示することもできます。dsp.SpectrumEstimator
オブジェクトの出力はスペクトル データです。このデータを取得して、その後の処理で使用できます。
spectrumAnalyzer を使用したパワー スペクトルの推定
信号のパワー スペクトルを表示するには、spectrumAnalyzer
System object™ を使用できます。入力信号のダイナミクスを変更して、その変更が信号のパワー スペクトルに与える影響をリアルタイムに確認することができます。
初期化
正弦波ソースを初期化して正弦波とスペクトル アナライザーを生成し、信号のパワー スペクトルを表示します。入力正弦波は 2 つの周波数 (1 つは 1000 Hz、もう 1 つは 5000 Hz) をもちます。2 つの dsp.SineWave
オブジェクトを作成します。1 つは 1000 Hz 正弦波を生成し、もう 1 つは 5000 Hz 正弦波を生成するオブジェクトです。
Fs = 44100; Sineobject1 = dsp.SineWave(SamplesPerFrame=1024,PhaseOffset=10,... SampleRate=Fs,Frequency=1000); Sineobject2 = dsp.SineWave(SamplesPerFrame=1024,... SampleRate=Fs,Frequency=5000); SA = spectrumAnalyzer(SampleRate=Fs,... SpectrumType="power",PlotAsTwoSidedSpectrum=false,... ChannelNames={'Power spectrum of the input'},YLimits=[-120 40],... ShowLegend=true);
スペクトル アナライザーはフィルター バンク アプローチを使用して、信号のパワー スペクトルを計算します。
推定
ストリーミングを行い、信号のパワー スペクトルを推定します。for
ループを作成して、5000 回の反復を実行します。各反復では、各正弦波の 1024 サンプル (1 フレーム) をストリーミングして、各フレームのパワー スペクトルを計算します。入力信号を生成するために、2 つの正弦波を追加します。結果として得られる信号は 2 つの周波数 (1 つは 1000 Hz、もう 1 つは 5000 Hz) をもつ正弦波です。ゼロ平均、標準偏差 0.001 のガウス ノイズを付加します。その後の処理で使用するスペクトル データを取得するには、オブジェクト関数 isNewDataReady
および getSpectrumData
を使用します。変数 data
には、スペクトルに関する追加の統計と共にスペクトル アナライザーに表示されるスペクトル データが含まれています。
data = []; for Iter = 1:1000 Sinewave1 = Sineobject1(); Sinewave2 = Sineobject2(); Input = Sinewave1 + Sinewave2; NoisyInput = Input + 0.001*randn(1024,1); SA(NoisyInput); if SA.isNewDataReady data = [data;getSpectrumData(SA)]; end end release(SA);
スペクトル アナライザーの出力では、2 つの明確なピーク (1 つは 1000 Hz、もう 1 つは 5000 Hz) を確認できます。
分解能帯域幅 (RBW) は、スペクトル アナライザーによって分解できる最小の周波数帯域幅です。既定では、spectrumAnalyzer
オブジェクトの RBWSource
プロパティは "auto"
に設定されます。このモードでは、RBW は周波数スパンと 1024 の比率です。両側スペクトルでは、この値は ですが、片側スペクトルでは です。この例のスペクトル アナライザーは片側スペクトルを示しています。そのため、RBW は (44100/2)/1024、つまり 21.53 Hz です。
のこの値を使用すると、1 つのスペクトルの更新を計算するのに必要な入力サンプル数 は、次の式で求められます。
この例では、 は 44100/21.53、つまり 2048 サンプルです。この値は、スペクトル アナライザーの下部にあるステータス バーで確認できます。
''自動'' モードで計算された は、適切な周波数分解能を提供します。
表示される 2 つの周波数を区別するには、2 つの周波数間の距離が少なくとも RBW でなければなりません。この例では、2 つのピーク間の距離は より大きい 4000 Hz です。そのため、ピークをはっきりと確認することができます。2 番目の正弦波の周波数を 1015 Hz に変更します。2 つの周波数間の差は より小さくなっています。
release(Sineobject2); Sineobject2.Frequency = 1015; for Iter = 1:1000 Sinewave1 = Sineobject1(); Sinewave2 = Sineobject2(); Input = Sinewave1 + Sinewave2; NoisyInput = Input + 0.001*randn(1024,1); SA(NoisyInput); end release(SA);
ピークは区別できません。
周波数分解能を向上させるには、 を 1 Hz まで減らします。
SA.RBWSource = "property"; SA.RBW = 1; for Iter = 1:3000 Sinewave1 = Sineobject1(); Sinewave2 = Sineobject2(); Input = Sinewave1 + Sinewave2; NoisyInput = Input + 0.001*randn(1024,1); SA(NoisyInput); end release(SA);
ズームすると、15 Hz 離れている 2 つのピークを区別できるようになりました。
周波数分解能を向上させると、時間分解能が低下します。周波数分解能と時間分解能の適切なバランスを維持するには、RBWSource
プロパティを "auto"
に変更します。
ストリーミング中に、入力のプロパティまたはスペクトル アナライザーのプロパティを変更して、スペクトル アナライザーの出力に与える影響を即座に確認することができます。たとえば、ループのインデックスが 1000 の倍数であるときに、2 番目の正弦波の周波数を変更します。
release(Sineobject2); SA.RBWSource = "auto"; for Iter = 1:5000 Sinewave1 = Sineobject1(); if (mod(Iter,1000) == 0) release(Sineobject2); Sineobject2.Frequency = Iter; Sinewave2 = Sineobject2(); else Sinewave2 = Sineobject2(); end Input = Sinewave1 + Sinewave2; NoisyInput = Input + 0.001*randn(1024,1); SA(NoisyInput); end release(SA);
ストリーミング ループの実行中、2 番目の正弦波のピークが反復値に従って変化することを確認できます。同様に、シミュレーションの実行中にスペクトル アナライザーの任意のプロパティを変更して、それに対応する出力の変化を確認できます。
単位間のパワーの変換
スペクトル アナライザーでは、パワー スペクトル密度を指定するための次の 3 つの単位が提供されています。[Watts/Hz]
、[dBm/Hz]
、および [dBW/Hz]
。対応するパワーの単位は [Watts]
、[dBm]
、および [dBW]
です。電気工学アプリケーションでは、信号の RMS を [Vrms]
または [dBV]
で表示することもできます。既定のスペクトル タイプは [dBm]
の [パワー] です。
ワット単位から dBW 単位および dBm 単位へのパワーの変換
[dBW]
単位のパワーは次で求められます。
[dBm]
単位のパワーは次で求められます。
振幅が 1 V の正弦波信号の場合、[Watts]
単位での片側スペクトルのパワーは次で求められます。
対応する dBm 単位のパワーは次で求められます。
ピークの検出でこの値を確認するには、スペクトル アナライザーのツールストリップの [測定値] タブで [ピークの検出] を有効にします。
ホワイト ノイズ信号の場合、スペクトルはすべての周波数でフラットになります。この例のスペクトル アナライザーは、[0 Fs/2] の範囲の片側スペクトルを表示します。分散 1e-4 のホワイト ノイズ信号の場合、単位帯域幅あたりのパワー (Punitbandwidth) は 1e-4 です。周波数範囲全体での [Watts]
単位のホワイト ノイズの合計パワーは次で求められます。
Fs は入力サンプル レートです。周波数ビンの数は、全帯域幅の RBW に対する比率です。片側スペクトルの場合、全帯域幅はサンプル レートの半分です。サンプル レートが 44100 Hz で RBW が 21.53 Hz の場合、[Watts]
単位のホワイト ノイズの合計パワーは 0.1024 W となります。
dBm 単位のホワイト ノイズのパワーは次を使用して計算できます。
ワット単位から dBFS 単位へのパワーの変換
スペクトル単位を [dBFS]
に設定し、フル スケール (FullScaleSource
) を "auto"
に設定すると、[dBFS]
単位のパワーは次のように計算されます。
ここでは以下のとおりです。
Pwatts はワット単位のパワー。
double および float の信号では、Full_Scale は入力信号の最大値。
固定小数点または整数の信号では、Full_Scale は表現できる最大の値。
手動のフル スケールを指定した場合 (FullScaleSource
を "property"
に設定)、[dBFS]
単位のパワーは次で求められます。
ここで、FS
は FullScale
プロパティで指定されたフル スケーリング係数です。
振幅が 1 V の正弦波信号の場合、[dBFS]
単位での片側スペクトルのパワーは次で求められます。
この値をスペクトル アナライザーで確認するには、次のコマンドを実行します。
Fs = 1000; % Sampling frequency sinef = dsp.SineWave(SampleRate=Fs,SamplesPerFrame=100); scope = spectrumAnalyzer(SampleRate=Fs,... SpectrumUnits="dBFS",PlotAsTwoSidedSpectrum=false) for ii = 1:100000 xsine = sinef(); scope(xsine) end
dBm 単位のパワーから Vrms 単位の RMS への変換
[dBm]
単位のパワーは次で求められます。
RMS の電圧は次で求められます。
前の例では、PdBm は 26.9897 dBm と同等です。Vrms は次のように計算されます。
0.7071 と等しくなります。
ピークの検出を使用してこの値を確認するには、次を行います。
スペクトル アナライザーのツールストリップの [アナライザー] タブで、[スペクトル]、[RMS] を選択します。
[測定値] タブで [ピークの検出] を有効にします。
dsp.SpectrumEstimator を使用したパワー スペクトルの推定
あるいは、dsp.SpectrumEstimator
System object を使用して信号のパワー スペクトルを計算することもできます。スペクトル推定器の出力を取得してデータを保存し、その後の処理で使用することができます。Estimation
ライブラリ内の他のオブジェクトを表示するには、MATLAB® コマンド プロンプトで help dsp
と入力して Estimation
をクリックします。
初期化
spectrumAnalyzer
の使用に関する前のセクションと同じソースを使って、パワー スペクトルを推定します。入力正弦波は 2 つの周波数 (1 つは 1000 Hz、もう 1 つは 5000 Hz) をもちます。dsp.SpectrumEstimator
を初期化し、フィルター バンク アプローチを使用して信号のパワー スペクトルを計算します。dsp.ArrayPlot
オブジェクトを使用して信号のパワー スペクトルを表示します。
Fs = 44100; Sineobject1 = dsp.SineWave(SamplesPerFrame=1024,PhaseOffset=10,... SampleRate=Fs,Frequency=1000); Sineobject2 = dsp.SineWave(SamplesPerFrame=1024,... SampleRate=Fs,Frequency=5000); SpecEst = dsp.SpectrumEstimator(Method='Filter bank',... PowerUnits='dBm',SampleRate=Fs,FrequencyRange='onesided'); ArrPlot = dsp.ArrayPlot(PlotType='Line',ChannelNames={'Power spectrum of the input'},... YLimits=[-80 30],XLabel='Number of samples per frame',YLabel='Power (dBm)',... Title='One-sided power spectrum with respect to samples');
推定
ストリーミングを行い、信号のパワー スペクトルを推定します。for
ループを作成して、5000 回の反復を実行します。各反復では、各正弦波の 1024 サンプル (1 フレーム) をストリーミングして、各フレームのパワー スペクトルを計算します。ゼロ平均、標準偏差 0.001 のガウス ノイズを入力信号に付加します。
for Iter = 1:5000 Sinewave1 = Sineobject1(); Sinewave2 = Sineobject2(); Input = Sinewave1 + Sinewave2; NoisyInput = Input + 0.001*randn(1024,1); PSoutput = SpecEst(NoisyInput); ArrPlot(PSoutput); end
フィルター バンク アプローチを使用すると、スペクトル推定の分解能は高くなり、ピークはスペクトル漏れのない正確なものになります。
x 軸を変換して周波数を表現
既定では、配列プロットはフレームあたりのサンプル数についてのパワー スペクトル データを表示します。x 軸の点の数は入力フレームの長さと等しくなります。スペクトル アナライザーは周波数に対してパワー スペクトル データをプロットします。片側スペクトルの場合、周波数は [0 Fs/2] の範囲で変化します。両側スペクトルの場合、周波数は [-Fs/2 Fs/2] の範囲で変化します。配列プロットの x 軸をサンプルベースから周波数ベースに変換するには、以下を実行します。
[Array Plot] ツールストリップの [スコープ] タブで、[設定] をクリックします。
片側スペクトルの場合: [データと座標軸] で、[サンプルのインクリメント] を に設定し、[X オフセット] を 0 に設定します。
両側スペクトルの場合: [データと座標軸] で、[サンプルのインクリメント] を に設定し、[X オフセット] を に設定します。
この例では、スペクトルは片側であるため、[サンプルのインクリメント] と [X オフセット] はそれぞれ 44100/1024 と 0 に設定されます。周波数を kHz 単位で指定するには、[サンプルのインクリメント] を 44.1/1024 に設定します。
ArrPlot.SampleIncrement = (Fs/1000)/1024; ArrPlot.XLabel = 'Frequency (kHz)'; ArrPlot.Title = 'One-sided power spectrum with respect to frequency'; for Iter = 1:5000 Sinewave1 = Sineobject1(); Sinewave2 = Sineobject2(); Input = Sinewave1 + Sinewave2; NoisyInput = Input + 0.001*randn(1024,1); PSoutput = SpecEst(NoisyInput); ArrPlot(PSoutput); end
ライブ処理
dsp.SpectrumEstimator
オブジェクトの出力にはスペクトル データが含まれており、その後の処理で使用できます。このデータはリアルタイムで処理することも、ワークスペースに格納することもできます。