Main Content

spa

スペクトル解析を使用した固定の周波数分解能による周波数応答の推定

説明

G = spa(data) は、時間領域または周波数領域のデータ data から、不確かさとノイズ スペクトルと共に周波数応答を推定します。data が時系列の場合、spa(data) は不確かさと共に出力パワー スペクトルを返します。spa は、ハン ウィンドウを使用して、0 (この値は除く) から π までの等間隔な 128 個の周波数値でスペクトルを計算します。

G = spa(data,winSize,freq) は、サイズ winSize のハン ウィンドウを使用して、freq に格納された周波数で周波数応答を推定します。

G = spa(data,winSize,freq,maxSize) は、各セグメントに含まれる要素が maxSize よりも少なくなるように入出力データをセグメントに分割します。この構文は計算性能を改善するために使用します。

すべて折りたたむ

iddata オブジェクト z3 の入出力データの周波数応答を推定します。既定の固定分解能である、0 (この値は除く) から π までの対数的に等間隔な 128 個の周波数値を使用します。

load iddata3 z3; 
g = spa(z3); 
bode(g)

Figure contains 2 axes objects. Axes object 1 with title From: u1 To: y1 contains an object of type line. This object represents g. Axes object 2 contains an object of type line. This object represents g.

対数的に等間隔なベクトル f を生成します。

f = logspace(-2,pi,128);

入出力データ z3 の周波数応答を推定します。ウィンドウ サイズを [] と指定して、既定のラグ ウィンドウ サイズを取得します。

load iddata3 z3;
g = spa(z3,[],f);

ボード線図と外乱スペクトルを標準偏差 3 の信頼区間でプロットします。

h = bodeplot(g);
showConfidence(h,3)

Figure contains 2 axes objects. Axes object 1 with title From: u1 To: y1 contains an object of type line. This object represents g. Axes object 2 contains an object of type line. This object represents g.

figure
h = spectrumplot(g);
showConfidence(h,3)

Figure contains an axes object. The axes object with title From: e@y1 To: y1 contains an object of type line. This object represents g.

入力引数

すべて折りたたむ

入出力データ。iddata オブジェクト、または複素数値を格納できる idfrd オブジェクトとして指定します。data には、出力のみをもつ時系列データも格納できます。

ハン ウィンドウのサイズ。ラグ サイズとも呼ばれます。スカラー整数として指定します。既定では、ウィンドウ サイズは min(length(data)/10,30) に設定されます。

スペクトル応答を推定する周波数。ラジアン/TimeUnit 単位の行ベクトルとして指定します。ここで、TimeUnitdataTimeUnit プロパティを表します。既定では、freq は範囲 (0,π] の対数的に等間隔な 128 個の値のベクトルに設定されます。離散時間モデルの場合は、ナイキスト周波数までの範囲内で freq を設定します。

data 内のセグメントの最大サイズ。正の整数として指定します。この引数を省略すると、関数はデータをセグメント化せずに、data 内のすべてのデータ セットを使用して推定を実行します。

出力引数

すべて折りたたむ

不確かさとノイズ スペクトルを含む周波数応答。idfrd オブジェクトとして指定されます。時系列データの場合、G は推定されたスペクトルと標準偏差になります。

推定結果と使用されたオプションに関する情報は、モデルの Report プロパティに格納されます。Report には次のフィールドがあります。

Report のフィールド説明
Status

モデルのステータスの概要。モデルが構築によって作成されたものか推定によって取得されたものかを示します。

Method

使用された推定コマンド。

windowSize

ハン ウィンドウのサイズ。

DataUsed

推定に使用されたデータの属性。次のフィールドをもつ構造体。

フィールド説明
Name

データ セットの名前。

Type

データ型。

Length

データ サンプルの数。

Ts

サンプル時間。これは、data.Ts と等価です。

InterSample

入力サンプル間動作。次の値のいずれかになります。

  • 'zoh' — ゼロ次ホールドにより、サンプル間で区分的に一定な入力信号を維持。

  • 'foh' — 1 次ホールドにより、サンプル間で区分的に線形な入力信号を維持。

  • 'bl' — 帯域幅を制限した動作により、ナイキスト周波数を超える連続時間入力信号のパワーがゼロになるように指定。

InterSample の値は離散時間モデルの推定結果には影響しません。

InputOffset

推定時に時間領域の入力データから削除されたオフセット。

OutputOffset

推定時に時間領域の出力データから削除されたオフセット。

詳細

すべて折りたたむ

周波数応答関数

"周波数応答関数" は、正弦波入力に対するシステムの定常状態応答を表します。線形システムの場合、特定の周波数の正弦波入力に対する結果は、同じ周波数をもち、異なる振幅と位相をもつ正弦波の出力になります。周波数応答関数は、振幅の変化と位相シフトを周波数の関数として表します。

周波数応答関数の理解を深めるために、次のように記述された線形動的システムについて考えてみます。

y(t)=G(q)u(t)+v(t)

ここで、u(t) と y(t) はそれぞれ入力信号と出力信号です。G(q) は、入力を出力に渡すシステム ダイナミクスを取得したもので、システムの伝達関数と呼ばれます。G(q)u(t) という表記は、次の演算を表しています。

G(q)u(t)=k=1g(k)u(tk)

q は "シフト演算子" で、次の方程式で定義されます。

G(q)=k=1g(k)qk      q1u(t)=u(t1)

G(q) は、単位円 G(q=e) で評価されるときの "周波数応答関数" です。

G(q=e) と出力ノイズ スペクトル Φ^v(ω) の組み合わせにより、周波数領域のシステムの記述が構成されます。

ブラックマン-テューキー法を使用して推定される周波数応答関数は、次の方程式で与えられます。

G^N(eiω)=Φ^yu(ω)Φ^u(ω)

この場合、^ はおおよその量を表します。この方程式の導出については、[1] の時間領域または周波数領域のノンパラメトリック手法に関する章を参照してください。

出力ノイズ スペクトル

出力ノイズ スペクトル (v(t) のスペクトル) は次の方程式で与えられます。

Φ^v(ω)=Φ^y(ω)|Φ^yu(ω)|2Φ^u(ω)

このノイズ スペクトルの方程式は、線形関係 y(t)=G(q)u(t)+v(t) が成り立ち、u(t) が v(t) から独立していて、スペクトル間に次の関係が成り立つと仮定して導出されたものです。

Φy(ω)=|G(eiω)|2Φu(ω)+Φv(ω)

Φyu(ω)=G(eiω)Φu(ω)

ここで、ノイズ スペクトルは次の方程式で与えられます。

Φv(ω)τ=Rv(τ)eiwτ

Φ^yu(ω) は出力と入力の相互スペクトル、Φ^u(ω) は入力のスペクトルです。

あるいは、外乱 v(t) をフィルター処理されたホワイト ノイズとして記述できます。

v(t)=H(q)e(t)

ここで、e(t) は分散 λ のホワイト ノイズであり、ノイズ パワースペクトルは次の方程式で与えられます。

Φv(ω)=λ|H(eiω)|2

アルゴリズム

spa は、次の手順に従ってブラックマン-テューキー スペクトル解析法を適用します。

  1. 共分散と相互共分散を u(t) と y(t) から計算します。

    R^y(τ)=1Nt=1Ny(t+τ)y(t)R^u(τ)=1Nt=1Nu(t+τ)u(t)R^yu(τ)=1Nt=1Ny(t+τ)u(t)

  2. 共分散と相互共分散のフーリエ変換を計算します。

    Φ^y(ω)=τ=MMR^y(τ)WM(τ)eiωτΦ^u(ω)=τ=MMR^u(τ)WM(τ)eiωτΦ^yu(ω)=τ=MMR^yu(τ)WM(τ)eiωτ

    ここで、WM(τ) は幅 (ラグ サイズ) M のハン ウィンドウです。M を指定して推定の周波数分解能を制御できます。2π/M ラジアン/サンプル時間とほぼ等しくなります。

    既定では、0 (この値は除く) から π までの等間隔な 128 個の周波数値がこの演算に使用されます。ここで、w = [1:128]/128*pi/Ts であり、Ts はそのデータ セットのサンプル時間です。ハン ウィンドウの既定のラグ サイズは M = min(length(data)/10,30) です。既定の周波数の場合は高速フーリエ変換 (FFT) が演算に使用され、ユーザー定義の周波数よりも効率が高くなります。

    メモ

    M = γ は、[1] の表 6.1 にあります。標準偏差は、[1] の 184 ページと 188 ページに記載されています。

  3. 周波数応答関数 G^N(eiω) と出力ノイズ スペクトル Φ^v(ω) を計算します。

    G^N(eiω)=Φ^yu(ω)Φ^u(ω)

    Φv(ω)τ=Rv(τ)eiwτ

spectrum は、出力と入力の両方のチャネルのスペクトル行列です。つまり、z = [data.OutputData, data.InputData] の場合、spectrum には z の行列値のパワー スペクトルがスペクトル データとして格納されます。

S=m=MMEz(t+m)z(t)WM(Ts)exp(iωm)

ここで、' は複素共役転置です。

参照

[1] Ljung, Lennart. System Identification: Theory for the User. 2nd ed. Prentice Hall Information and System Sciences Series. Upper Saddle River, NJ: Prentice Hall PTR, 1999.

バージョン履歴

R2006a より前に導入