Main Content

spa

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

説明

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 (この値は除く) から π までの対数的に等間隔な 128 個の周波数値を使用します。

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

Figure contains 2 axes objects. Axes object 1 with title From: u To: y, ylabel Magnitude (dB) contains an object of type line. This object represents g. Axes object 2 with ylabel Phase (deg) contains an object of type line. This object represents 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 contains 2 axes objects. Axes object 1 with title From: u1 To: y1, ylabel Magnitude (dB) contains an object of type line. This object represents g. Axes object 2 with ylabel Phase (deg) 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, ylabel Power (dB) contains an object of type line. This object represents g.

入力引数

すべて折りたたむ

等間隔にサンプリングされた推定データ。以降のセクションで説明しているとおり、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 — Ns 行 Nu 列の行列。Nu は入力の数です。

  • y — Ns 行 Ny 列の行列。Ny は出力の数です。

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

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

データ オブジェクト

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

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

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

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

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

名前と値の引数

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

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

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

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

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

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

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

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

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

出力引数

すべて折りたたむ

不確かさとノイズ スペクトルを含む周波数応答。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 より前に導入

すべて展開する