離散化された入出力デ​ータから、周波数応答​を得る方法

離散化された入出力データから、周波数応答を得ようとしています。

まずはシミュレーションにより得たデータによる周波数応答と、シミュレーションにおける伝達関数のボード線図を比較して一致するかどうかを確かめようとしているのですが、固有角周波数が一致していません。以下、解析対象と解析手順について説明します。

■解析対象

シミュレーションにおける制御対象の運動方程式：m* d^2 r/(dt)^2 =F

シミュレーションにおいて使用する制御則（微分先行型PD制御）：F = kp*(rd - r) - kd*dr/dt

○以下パラメータ

・質点の質量：m=1 [kg]

・比例ゲイン：kp = 100 [N/m]

・微分ゲイン：kd = 25 [N/m]

・固有角周波数：ωｎ = 10 [rad/s]

・減衰係数：ζ = 1

・サンプル周波数：Fs = 1000 [Hz]

○以上のパラメータの場合、ボード線図は図1のように表されます。

図1．微分先行型PD制御を適応した場合のボード線図

■周波数応答解析の手順

上記のパラメータの解析対象において、以下の入出力データを用いて周波数応答の解析を試みました。

・入力…目標値：rd =1 [m]

・出力…質点の位置：r [m]

ただし、入力であるrdはステップ時間を1 [s]としたステップ入力として与えました。

使用した関数は[Txy,f] = tfestimate(x,y,window,noverlap,nfs,Fs)です。結果を図2に示します。

実際のプログラムコードを以下に記します。

Fs=1000%サンプル周波数
nfs = 16384%サンプリング点数（分解能 = 約0.06）
N = 10%窓長をnfsの何分の1にするか
wind = hann(nfs/N)%ハニング窓（窓の長さ調整）
overlap = nfs/(2*N)%オーバラップ値。窓長の半分までオーバラップする
[Txy,f] = tfestimate(rd_PD_simu,r_PD_simu,wind,overlap,nfs,Fs);

semilogx(f,20*log10(abs(Txy)))%対数グラフのプロット

xlabel('ω [rad/s]')
ylabel('20*log10|A|')
hold on

%% 微分先行型PDの伝達関数をボード線図にプロット
Gs_PD = tf([100],[1 20 100]);
bode(Gs_PD,'r')
set(gca, 'FontSize', 13)

以上の手順を踏み、結果は図2のようになりました。

図2.　解析結果と伝達関数の比較

図2において、解析結果である青のグラフが赤のグラフと一致する必要がありますが、一致できていません。（ただし、10＾2～10＾3の範囲で値が急激に減少しているのはサンプル周波数Fsの影響です。）固有角周波数が青のグラフにおいては10＾1となっていないことが不自然です。なぜこのような結果になってしまっているのか、わかる方はおられるでしょうか？

Answer 1

[Naoya]

tfestimate関数の戻り値の f (周波数) の単位は、 Hz になります。

一方で、 tf 関数からの連続システムにおける ボード線図の周波数の単位は、そのボード線図の応答から判断すると、rad/sec のようです。

恐らく、周波数の単位の不一致分のシフト量の誤差と推測します。

例えば、rad/sec にあわせるのであれば、

semilogx(f,20*log10(abs(Txy)))

を

semilogx(f/(2*pi),20*log10(abs(Txy))) 

にして確認していただけますでしょうか？

[Naoya]

すみません、上の件ですが、 Hz ⇒ rad/sec ですので、

semilogx(f*(2*pi),20*log10(abs(Txy))) 

となります。

[Takahiro OHASHI]

無事解決いたしました。

戻り値fの単位を誤認していたんですね。

ありがとうございます！