データ取得

まずは適切な入力を作成して入出力のデータを得る。

Chirp信号作成

簡単に作れて，かつそこそこの同定をするならChirp信号を用いる。
simulinkのブロックはなんか挙動が思ったように行かなかったので，Mfunctionで作成。

Linearな信号で良ければ1行で作成可能。
下記wikipediaの記事を参照。
https://ja.wikipedia.org/wiki/%E3%83%81%E3%83%A3%E3%83%BC%E3%83%97%E4%BF%A1%E5%8F%B7

u: clockから得た時間
f0: 最小の周波数
f1: 最大の周波数
Tchirp: 一度のsweep時間
phi0: 初期位相（0でよさそう）


sin(phi0+2*pi*(f0*rem(u,Tchirp)+(f1-f0)/2/Tchirp*rem(u,Tchirp)^2))

こんなかんじ。rem関数を使って周期的な動作を可能にしているのがミソ。

解析

そんでもって得たデータをin,outとする。

以下のがその一例。
(https://qiita-image-store.s3.amazonaws.com/0/188863/1243447c-53d8-b6a2-4eee-6a1d2d42dfa3.png)

周波数応答の同定

freqは当然周波数だが，winsizeは得られたデータをどのサイズで窓関数に打ち込むかを与える。
感覚的にはTchirp*freqで良さそうだが，間違ってたら指摘してください。


[PXX,FREQ] = tfestimate(in,out,winsize,[],[],freq);
Txx = frd(PXX,FREQ*2*pi,1/freq);

2行目でFREQに2*piをかけているのはHzからradに変換する必要があるため。どこかパラメータ設定すればこんなことしなくて済んだ気もする。

ボード線図を引いて見ればキレイに表示されていることと思う。


Bode(Txx)

また，Coherenceを取得するのもコマンド一つで済む。

[COHER,cFREQ] = mscohere(in,out,winsize,[],[],freq);
semilogx(cFREQ,COHER);

coherenceが小さいほどあんまり信用できないデータということになる。

お手軽フィッティング＆システム同定

invfreqsという関数を用いることで簡単にフィッティングができる。学生はmatlabに課金すべき。

コマンドの使い方は以下の通りで先程tfestimateで生成したPXXとFREQを用いる。（必要に応じてradに変換すること！）

[num den] = invfreqs(PXX,FREQ*2*pi,nd,dd,weight,iter)

3,4番目のnd,ddは分子と分母の字数で1次/2次なら,1,2と綴る。![octave-online-line-30.png]
weightは重み関数である。Coherenceから作成するとスムース。
iterはイテレーション回数。大体10回やっておけばまぁまぁ近い値になる。

重み関数作り方例

周波数帯がminf～maxfの間の情報のみかつ，Coherenceが0.95以上の周波数応答を用いてfittingするとする。
こんなかんじに書ける。簡単でしょう。

weight = zeros(size(COHER));
weight(COHER>0.95) = 1;
weight(cFREQ >= fmax)= 0;
weight(cFREQ <= fmin)= 0;

Sample Code（バグあり）

今matlabが消えたPCでやってるので後でちゃんと内容を確認する。あとで。

%% simulation set
close all
freq = 1000; %1000hz

%% 0. Define system
% system (ground truth) 
sys = tf([0.6],[0.02,10]);
figure;bode(sys);

%% 1. read data and set in/out data example

% make chirp signals
time = 0:1/freq:10;
CHP=chirp(time,0.5,2,50);
out=lsim(sys,CHP,time);
in=CHP;

%%load('result.data');
%time = result(:,1);
%in = result(:,2);
%out = result(:,3);

figure(1); plot(time,in,time,out);
xlabel('time[s]')
ylabel('in[A]/out[rad/s]')
legend('In','Out')
grid on;

%% 2. Get bode plot with tfestimate
windowsize = 2000;% this is not nessesary
[PXX,FREQ] = tfestimate(in,out,windowsize,[],[],freq);
Txx = frd(PXX,FREQ,1/freq);

figure(2); 
%opt = bodeoptions('cstprefs');
%opt.PhaseWrapping = 'on';
%opt.PhaseWrappingBranch = -180;
%bode(Txx,opt);
bode(Txx);
title('Freq responce')

%% 3. Get coherence mscoherence
[COHER,cFREQ] = mscohere(in,out,windowsize,[],[],freq);

figure(3); semilogx(cFREQ,COHER);
title('coherent')
xlabel('Freq [Hz]')
%% 4. fitting1
% Chirp Sweep Range
fmax = 50;
fmin = 0.5;

% choose weight from coherence
weight = zeros(size(COHER));
weight(COHER>0.95) = 1;
% use only 0.5--50 Hz
weight(cFREQ >= fmax)= 0;
weight(cFREQ <= fmin)= 0;

% estimate transfer function !
[num den]=invfreqs(PXX,2*pi*FREQ,0,1,weight,10);

% figure;bode(tf(num,den));
fitsys = tf(num,den);

figure(4);bode(sys,fitsys);
legend('ground truth','Estimated')
xlim([fmin fmax])

MATLABで簡単システム同定（簡易）