■実現したいこと
今回解決したいこととは別に、自身が実現したいことを最初に記しておきます。
私がしたいことは、Matlab Simulinkでのシミュレーションにおいて、0.10～200 [Hz]の周波数応答（入出力から得た伝達関数の推定）を見るということです。
■現状問題となっていること
現状問題となっていることは、「解析が可能な範囲において、周波数応答が（低周波帯域において）正確に解析が出来ていない」ということです。
今回設定したパラメータなら解析が出来ると認識しているのですが、想定した周波数帯域まで正確に解析できていません。
■解析対象
今回の解析対象を説明します。
今回の解析対象は、Zero-Poleブロックで構成された、2次のバターワースフィルタを含む1次の微分器です。想定する微分器をラプラス変換したブロックを以下に示します。ただし、入力をx、出力をxの微分とした場合の図を表しています。
図1. 想定する微分器の構成
図1を模すために、Zero-Poleブロックを作成する上で必要なパラメータは以下のように設定しています。
・零点 ： 0
・極 ： [sqrt(1/2)*(-1 -i)*DiffLPF(1).omegac,sqrt(1/2)*(-1 +i)*DiffLPF(1).omegac]
・ゲイン ： (2*pi*omega_c)^n ※omega_c = 20 [Hz]（カットオフ周波数）、n = 2
(・絶対許容誤差：auto)
■周波数応答の解析方法とパラメータ設定
現状は周波数応答を得るために、"tfestimate"という関数を用いています。実際には次のように表される関数です。
[Txy,f] = tfestimate(x,y,wind,overlap,nfs,Fs)
これはFFTを用いて周波数応答を得ています。カッコ内のパラメータやその他のパラメータはそれぞれ、以下のように定めています。
・ｘ ：(Matlab Simulinkにおける)入力。
・ｙ ：(Matlab Simulinkにおける)出力。図におけるxの微分に相当。
・nfs = 131072 ：サンプリング点数。2^nで表す必要があります。
・Ｎ = 8 ：窓長の長さをサンプリング点数の何分の1にするかを決めます。
・wind = hann(nfs/N) ：何の窓関数をしようするか。現状はハニング窓を使用しています。カッコ内は窓長の長さを示しています。
・overlap = nfs/(2*N) ：オーバーラップ値。現状は窓長の半分の長さにしています。
・Fs = 1000 [Hz] ：サンプリング周波数。今回であればシミュレーションのものと一致させます。
・今回使用している入力と出力のデータ数 = 30000
■設定したパラメータから算出される解析範囲
以上の設定したパラメータより、得られる最大の周波数と最小の周波数は次のようになります。
ただし、
https://www.onosokki.co.jp/HP-WK/c_support/faq/fft_common/CF_9.htm
を参考にしました。
・最大の周波数 = Fs/2.56 = 1000/2.56 = 約390 [Hz]
・最小の周波数 = （最大の周波数*2.56*N）/nfs = （サンプリング周波数）/（窓長） = 1000*8/131072 = 約0.061 [Hz]
これは当初指定していた周波数の範囲を含みます。
■解析結果
以上の解析対象を解析した結果を添付した画像に示します。
ただし、青の線は解析を行った結果を表しており、赤の線は図1の想定した微分器の伝達関数をプロットしたものになります。
添付した画像から分かることは以下の2つです。
①高周波帯域においては1000 [rad/s]以降で誤差が生じている
・・・解析範囲は（1000/2.56 [Hz]）*2*pi = 約2500 [rad/s]までが最大
②低周波帯域においては10 [rad/s]以下で誤差が生じている
・・・解析範囲は（1000*8/131072 [Hz]）*2*pi = 約0.38 [rad/s]までが最小
■結果に対する自身の理解
解析結果に対して私がどのように解釈しているかを説明します。
①については、今回は微分器を対象としているため高周波帯域において誤差が出るのは理解できます。
ただ②については、何故今回のような結果になるのかわかりません。今回設定したパラメータにおいては約0.38 [rad/s]まで正確に測定できるはずですが、結果を見る限りは出来ていません。
以上、分かる方がいれば、教えていただければ幸いです。
