spa

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

構文

G = spa(data)

G = spa(data,winSize,freq)

G = spa(data,winSize,freq,maxSize)

sys = spa(___,Name,Value)

説明

G = spa(data) は、時間領域または周波数領域のデータ data から、不確かさとノイズ スペクトルと共に周波数応答を推定します。

data には、時間領域または周波数領域の入出力データまたは時系列データが含まれます。data は、timetable、数値行列のコンマ区切りのペア、iddata オブジェクト、または複素数値を含むことができる idfrd オブジェクトの形式にすることができます。

data が時系列の場合、つまり入力が含まれない場合、spa(data) は不確かさと共に出力パワー スペクトルを返します。spa は、ハン ウィンドウを使用して、0 (この値は除く) から π までの等間隔な 128 個の周波数値でスペクトルを計算します。

data が timetable の形式の場合、最後の変数が単一の出力変数であるものと解釈されます。この解釈を変更するには、名前と値の引数 InputName および OutputName を使用します。

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

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

sys = spa(___,Name,Value) は、1 つ以上の名前と値の引数で指定された追加のモデル オプションを使用します。

たとえば、sys = spa(data,'InputName',["u1","u3"],'OutputName',["y1","y4"]) を使用して入出力の信号変数名を指定します。

この構文では、前述の任意の入力引数の組み合わせで使用できます。

例

周波数応答の推定

timetable tt3 の入出力データの周波数応答を推定します。既定の固定分解能である、0 (この値は除く) から $$\pi$$ までの対数的に等間隔な 128 個の周波数値を使用します。

load sdata3 tt3; 
g = spa(tt3); 
bode(g)

指定した周波数での周波数応答の推定

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

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

行列 umat3 および ymat3 内の入出力データに対する周波数応答を推定します。ウィンドウ サイズを [] と指定して、既定のラグ ウィンドウ サイズを取得します。

load sdata3 umat3 ymat3;
g = spa(umat3,ymat3,[],f);

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

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

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

入力引数

data — 推定データ

timetable | 数値行列のペア | iddata オブジェクト | idfrd オブジェクト | timetable の cell 配列 | 行列の cell 配列

等間隔にサンプリングされた推定データ。以降のセクションで説明しているとおり、timetable、コンマ区切りの行列のペア、またはデータ オブジェクトとして指定します。

timetable

一定間隔の時間ベクトルを使用する timetable として data を指定します。data には入力チャネルと出力チャネルを表す変数が含まれます。

推定に timetable を使用する場合、すべての変数を使用するか、使用するチャネルのサブセットを指定できます。

複数実験データの場合は、データを timetable の Ne 行 1 列の cell 配列として指定します。Ne は実験数です。すべての実験のサンプル時間が一致していなければなりません。

推定に使用する個別の入出力チャネルを選択するには、名前と値の引数 InputName および OutputName を使用します。

時系列データの場合、出力変数のみが含まれている単一変数の timetable として data を指定するか、名前と値の引数 InputName および OutputName を使用して推定用の出力変数のみを指定します。

たとえば、sys = spa(data,'OutputName',"y2",'InputName',[]) は、出力変数 y2 を使用し、入力変数なしで、timetable data から時系列モデル sys を推定します。

コンマ区切りの行列ペア

等間隔にサンプリングされた入出力の時間領域の信号値が含まれた行列 u,y のコンマ区切りのペアとして data を指定します。行列ベースのデータはサンプル時間情報を提供しません。サンプル時間は 1 秒と仮定されます。連続時間システムでの行列ベースのデータの使用は推奨されません。

SISO システムの場合は、長さが Ns の列ベクトルとして u,y を指定します。ここで、Ns はサンプル数です。

MIMO システムの場合は、以下の次元をもつ行列ペアとして u,y を指定します。

• u — N_s 行 N_u 列の行列。N_u は入力の数です。

• y — N_s 行 N_y 列の行列。N_y は出力の数です。

複数実験データの場合は、u,y を 1 行 N_e 列の cell 配列のペアとして指定します。N_e は実験数です。

時系列システムの場合は、空の u、つまり [],y を指定します。

データ オブジェクト

推定データ オブジェクト。時間領域または周波数領域の iddata オブジェクト、または等間隔にサンプリングされた入出力値を含む idfrd オブジェクトとして指定します。データ オブジェクトは 1 つ以上の出力チャネルとゼロ個以上の入力チャネルをもつことができます。既定では、モデルのサンプル時間は推定データのサンプル時間に設定されます。

推定データ型の取り扱いの詳細については、Data Domains and Data Types in System Identification Toolboxを参照してください。

winSize — ウィンドウ サイズ

[] (既定値) | スカラー整数

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

freq — 周波数

行ベクトル

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

maxSize — 最大セグメント サイズ

250e3 (既定値) | 正の整数

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

名前と値の引数

引数のオプションのペアを Name1=Value1,...,NameN=ValueN として指定します。Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後に表示されなければなりませんが、ペアの順序は重要ではありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: sys = spa(data,'InputName',"u2")

InputName — 入力チャネル名

string | 文字ベクトル | string 配列 | 文字ベクトルの cell 配列

入力チャネル名。string、文字ベクトル、string 配列、または文字ベクトルの cell 配列として指定します。

データ ソースに timetable を使用している場合は、InputName の名前は timetable 変数のサブセットでなければなりません。

例: sys = spa(tt,__,'InputName',["u1" "u2"]) は、推定で使用する timetable tt からの入力チャネルとして変数 u1 および u2 を選択します。

OutputName — 出力チャネル名

string | 文字ベクトル | string 配列 | 文字ベクトルの cell 配列

出力チャネル名。string、文字ベクトル、string 配列、または文字ベクトルの cell 配列として指定します。

データ ソースに timetable を使用している場合は、OutputName の名前は timetable 変数のサブセットでなければなりません。

例: sys = spa(tt,__,'OutputName',["y1" "y3"]) は、推定で使用する timetable tt からの出力チャネルとして変数 y1 および y3 を選択します。

出力引数

G — 周波数応答とノイズ スペクトル

idfrd オブジェクト

不確かさとノイズ スペクトルを含む周波数応答。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)={\displaystyle \sum _{k=1}^{\infty }g(k)u(t-k)}$$

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

$$G(q)={\displaystyle \sum _{k=1}^{\infty }g(k)q^{-k}}{q}^{-1}u(t)=u(t-1)$$

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

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

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

$${\widehat{G}}_N\left(e^{i\omega }\right)=\frac{{\widehat{\Phi }}_{yu}\left(\omega \right)}{{\widehat{\Phi }}_u\left(\omega \right)}$$

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

出力ノイズ スペクトル

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

$${\widehat{\Phi }}_v\left(\omega \right)={\widehat{\Phi }}_y\left(\omega \right)-\frac{{\left|{\widehat{\Phi }}_{yu}\left(\omega \right)\right|}^2}{{\widehat{\Phi }}_u\left(\omega \right)}$$

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

$${\Phi }_y(\omega )={\left|G\left(e^{i\omega }\right)\right|}^2{\Phi }_u(\omega )+{\Phi }_v(\omega )$$

$${\Phi }_{yu}(\omega )=G\left(e^{i\omega }\right){\Phi }_u(\omega )$$

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

$${\Phi }_v(\omega )\equiv {\displaystyle \sum _{\tau =-\infty }^{\infty }{R}_v}(\tau )e^{-iw\tau }$$

$${\widehat{\Phi }}_{yu}(\omega )$$ は出力と入力の相互スペクトル、$${\widehat{\Phi }}_u(\omega )$$ は入力のスペクトルです。

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

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

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

$${\Phi }_v(\omega )=\lambda {\left|H\left(e^{i\omega }\right)\right|}^2$$

アルゴリズム

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

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

   $$\begin{array}{l}{\widehat{R}}_y\left(\tau \right)=\frac{1}{N}{\displaystyle \sum _{t=1}^{N}y\left(t+\tau \right)y\left(t\right)}\\ {\widehat{R}}_u\left(\tau \right)=\frac{1}{N}{\displaystyle \sum _{t=1}^{N}u\left(t+\tau \right)u\left(t\right)}\\ {\widehat{R}}_{yu}\left(\tau \right)=\frac{1}{N}{\displaystyle \sum _{t=1}^{N}y\left(t+\tau \right)u\left(t\right)}\end{array}$$

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

   $$\begin{array}{l}{\widehat{\Phi }}_y(\omega )={\displaystyle \sum _{\tau =-M}^M{\widehat{R}}_y}(\tau )W_M(\tau )e^{-i\omega \tau }\\ {\widehat{\Phi }}_u(\omega )={\displaystyle \sum _{\tau =-M}^M{\widehat{R}}_u}(\tau )W_M(\tau )e^{-i\omega \tau }\\ {\widehat{\Phi }}_{yu}(\omega )={\displaystyle \sum _{\tau =-M}^M{\widehat{R}}_{yu}}(\tau )W_M(\tau )e^{-i\omega \tau }\end{array}$$

   ここで、$$W_M(\tau )$$ は幅 (ラグ サイズ) 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. 周波数応答関数 $${\widehat{G}}_N\left(e^{i\omega }\right)$$ と出力ノイズ スペクトル $${\widehat{\Phi }}_v(\omega )$$ を計算します。

   $${\widehat{G}}_N\left(e^{i\omega }\right)=\frac{{\widehat{\Phi }}_{yu}\left(\omega \right)}{{\widehat{\Phi }}_u\left(\omega \right)}$$

   $${\Phi }_v(\omega )\equiv {\displaystyle \sum _{\tau =-\infty }^{\infty }{R}_v}(\tau )e^{-iw\tau }$$

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

$$S={\displaystyle \sum _{m=-M}^{M}Ez\left(t+m\right)z}{\left(t\right)}^{\prime }{W}_M\left(T_s\right)\exp \left(-i\omega m\right)$$

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

参照

[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 より前に導入

R2022b: 時間領域推定データが timetable 形式および行列形式で受け入れられる

ほとんどの推定関数、検証関数、解析関数、およびユーティリティ関数が、時間領域入出力データを、入力と出力両方のデータを含む単一の timetable の形式、または、入力データと出力データを別々に含む行列のペアの形式で受け入れるようになりました。これらの関数は、時間領域データと周波数領域データの両方について、iddata オブジェクトを引き続きデータ ソースとして受け入れます。

参考

etfe | freqresp | idfrd | spafdr | bode | spectrum

トピック

• What is a Frequency-Response Model?
• Estimate Frequency-Response Models at the Command Line
• Selecting the Method for Computing Spectral Models