ボード線図が分かって​いる場合の伝達関数の​求め方を教えてくださ​い

37 views (last 30 days)
Jiro Yamada
Jiro Yamada on 5 Apr 2019
Commented: Jiro Yamada on 24 Apr 2019
以下のコードで得られるボード線図の伝達関数を求めたいです。
E=7.06*10^10;
B=3.0*10^(-2);
H=1.5*10^(-3);
l=1.0;
I=(B*H^3)/12;
p=2680;
A=B*H;
xL=0.05;
w=0.1:1/1000:100;
wr=w*2*pi;
k=((p*A.*wr.^2)./(E*I)).^(1/4);
c=exp(1i.*k.*xL);
F1=c;
ReF1=real(F1);
ImF1=imag(F1);
rF1=sqrt(ReF1.^2+ImF1.^2);
rF1dB=20*log10(rF1);
sitaF1=(atan(ImF1./ReF1))*(180/pi);
figure(1)
subplot(2,1,1)
semilogx(w,rF1dB,'r')
ylim([-3 5])
ylabel('Gain [dB]')
semilogx(w,sitaF1,'r')
xlabel('Frequency [Hz]')
ylabel('Phase [degree]')
現在はinvfreqsを用いた以下のコードを使って伝達関数を求めているのですが、希望のボード線図となるようなものは得られていません。
[F1b,F1a]=invfreqs(F1,wr,1,0);
sys_s_F1=tf(F1b,F1a);
どのようにすればfigure(1)のボード線図となる伝達関数を求めることができるのでしょうか?

Accepted Answer

Hiroumi Mita
Hiroumi Mita on 8 Apr 2019
掲題の問題について
以下のようにシステム同定の手法でモデルを推定しますと
一応、ほぼにたようなボード線図を書くことはできます。
しかし、得られた伝達関数は複素係数のものになり
現実性は乏しく、無理やり数値計算で合わせこんだ印象があります。
もともと推定したいモデルのボード線図は
周波数帯域において、ゲイン特性が一定数(0dB)、位相は
周波数が上がるとじりじり上がるもので、現実的に
このようなモデルがあるのかどうか、必要性はあるのか、
数式表現できるものなのか大いに疑問があります。
この問題の前提条件を再確認されることをお勧めします。
system_ident.bmp
E=7.06*10^10;
B=3.0*10^(-2);
H=1.5*10^(-3);
l=1.0;
I=(B*H^3)/12;
p=2680;
A=B*H;
xL=0.05;
w=0.1:1/1000:100;
wr=w*2*pi;
k=((p*A.*wr.^2)./(E*I)).^(1/4);
c=exp(1i.*k.*xL);
F1=c;
ReF1=real(F1);
ImF1=imag(F1);
rF1=sqrt(ReF1.^2+ImF1.^2);
rF1dB=20*log10(rF1);
sitaF1=(atan(ImF1./ReF1))*(180/pi);
figure(1)
subplot(2,1,1)
semilogx(w,rF1dB,'r')
ylim([-3 5])
ylabel('Gain [dB]')
subplot(2,1,2)%
semilogx(w,sitaF1,'r')
xlabel('Frequency [Hz]')
ylabel('Phase [degree]')
%システム同定
[F1b,F1a]=invfreqs(F1,wr,1,0);
sys_s_F1=tf(F1b,F1a);
mag=rF1
phase = (atan(ImF1./ReF1))*(180/pi);
response = mag.*exp(1j*phase*pi/180);
fr_data = idfrd(response,wr,1/(100*2));
% Import mydata
% State space model estimation
Options = n4sidOptions;
Options.Display = 'on';
Options.N4Weight = 'CVA';
Options.N4Horizon = [15 15 15];
ss1 = n4sid(fr_data, 9, 'Ts', 0, Options)
tf(ss1)%伝達関数表示
figure;bode(ss1,{0,100*2*pi})%ボード線図
  2 Comments
Jiro Yamada
Jiro Yamada on 17 Apr 2019
度々すみません。
位相が下がっていくこちらならどうでしょうか。
同じコードで試したのですが、なかなかにかけ離れたものが得られてしまいました。
どこを調節すればよいかなども教えていただけると幸いです。
たくさんお願いして申し訳ありませんが、よろしくお願いします。
%----------パラメータ設定-----------%
E=7.06*10^10; %縦弾性係数[N/m^2]
B=3.0*10^(-2); %はりの幅[m]
H=1.5*10^(-3); %はりの厚さ[m]
l=1.0; %はりの全長[m]
I=(B*H^3)/12; %断面二次モーメント
p=2680; %密度[kg/m^3]
A=B*H; %はりの断面積
xL=0.05; %センサ間隔[m]
w=0.1:1/1000:100; %[Hz]
wr=w*2*pi; %[rad/s]
k=((p*A.*wr.^2)./(E*I)).^(1/4);
L10=0.3;
L20=0.7;
L30=1.0;
L21=L20-L10;
L31=L30-L10;
L32=L30-L20;
Ls1=0.1;
%----各パラメータ--------------------
delta=ijTxy(3,3,k,L30).*ijTxy(4,4,k,L30)-ijTxy(3,4,k,L30).*ijTxy(4,3,k,L30);
alpha11=ijTxy(4,4,k,L30).*ijTxy(3,3,k,L31)-ijTxy(3,4,k,L30).*ijTxy(4,3,k,L31);
alpha12=ijTxy(4,4,k,L30).*ijTxy(3,3,k,L32)-ijTxy(3,4,k,L30).*ijTxy(4,3,k,L32);
alpha21=ijTxy(3,3,k,L30).*ijTxy(4,3,k,L31)-ijTxy(4,3,k,L30).*ijTxy(3,3,k,L31);
alpha22=ijTxy(3,3,k,L30).*ijTxy(4,3,k,L32)-ijTxy(4,3,k,L30).*ijTxy(3,3,k,L32);
beta1=ijTxy(1,3,k,Ls1)/4+1j*ijTxy(2,3,k,Ls1)./(4*k)-ijTxy(3,3,k,Ls1)./(4*k.^2)-1j*ijTxy(4,3,k,Ls1)./(4*k.^3);
beta2=ijTxy(1,4,k,Ls1)/4+1j*ijTxy(2,4,k,Ls1)./(4*k)-ijTxy(3,4,k,Ls1)./(4*k.^2)-1j*ijTxy(4,4,k,Ls1)./(4*k.^3);
gamma1=k.^3.*ijTxy(1,3,k,L10)+1j.*k.^2.*ijTxy(2,3,k,L10)-k.*ijTxy(3,3,k,L10)-1j.*ijTxy(4,3,k,L10);
gamma2=k.^3.*ijTxy(1,4,k,L10)+1j.*k.^2.*ijTxy(2,4,k,L10)-k.*ijTxy(3,4,k,L10)-1j.*ijTxy(4,4,k,L10);
GN=delta.*(alpha12.*gamma1+alpha22.*gamma2);
GD=(alpha11.*alpha22-alpha12.*alpha21).*(beta2.*gamma1-beta1.*gamma2)+delta.*k.*(alpha12.*beta1+alpha22.*beta2);
G=GN./GD;
%希望の周波数特性
ReG=real(G);
ImG=imag(G);
rG=sqrt(ReG.^2+ImG.^2); %[dB]
rGdB=20*log10(rG);
sitaG=unwrap(angle(G))*(180/pi); %[degree]
figure(1)  %希望の周波数特性
subplot(2,1,1)
semilogx(w,rGdB,'k')
xlim([0.1 100])
ylabel('Gain [dB]')
subplot(2,1,2)
semilogx(w,sitaG,'k')
xlim([0.1 100])
xlabel('Frequency [Hz]')
ylabel('Phase [degree]')
%ここから教えていただいたところ
mag=rG;
phase = (atan(ImG./ReG))*(180/pi);
response = mag.*exp(1j*phase*pi/180);
fr_data = idfrd(response,wr,1/(100*2));
% Import mydata
% State space model estimation
Options = n4sidOptions;
Options.Display = 'on';
Options.N4Weight = 'CVA';
Options.N4Horizon = [15 15 15];
ss1 = n4sid(fr_data, 9, 'Ts', 0, Options);
tf(ss1)%伝達関数表示
figure(2)
bode(ss1,{0,100*2*pi})%ボード線図

Sign in to comment.

More Answers (2)

Hiroumi Mita
Hiroumi Mita on 17 Apr 2019
ボード線図において、例えば次のような特徴を持つものは
無駄時間が該当します。
#1. ゲイン特性が周波数帯域において一様0[dB]
#2. 位相特性は周波数が高くなると遅れる
このサンプルが無駄時間1[s]のボード線図です。
s=tf('s')
sys=exp(-s)
figure;bode(sys)
無駄時間は入力が入っても一定時間(無駄時間)無反応なものを
言います。
伝達関数は
sys=exp(-T*s)
T:無駄時間
です。
もし、このシステムを対象にするなら、出力が何秒遅れて入力に反応したかを
調べればよく、ボード線図から推定するのは意味がありません。
何らかの動特性を持つシステムは、無駄時間のみ表現されるものはないと思われます。
動特性(時間による変化)+無駄時間が一般的です。
基本的に、ボード線図からモデルを推定するとき、周波数の変化に対する
ゲイン特性の変化、位相の変化を必要とします。
そこで、示されたボード線図は、何らかの実験の不手際により、
ゲイン特性が正しく取れなかったか?、元の仮定では動特性を無視したのか?
そんなことが考えられるわけです。
動特性は示された系に無いのか?もし、本当にないのならそれはボード線図で扱う領域ではないかもしれません。
そのあたりを再考いただければと思います。
  1 Comment
Jiro Yamada
Jiro Yamada on 18 Apr 2019
解説までしていただきありがとうございます。
Simulinkで変位センサの信号を通すデジタルフィルタを作りたくて伝達関数の推定を行っております。
希望のボード線図は実験で得られたものではなく、数値計算を基に求めたものになります(最初のコードの13行目にあるF1=cから)。
教えていただいたことを参考に考えてみます。

Sign in to comment.


Shoumei
Shoumei on 24 Apr 2019
位相が右上がりではなく右下がりで良ければこんな特性のフィルタはありますが。
ギターのエフェクトでよく使われるフェイザーのオールパスフィルタです。
APF.jpg
  1 Comment
Jiro Yamada
Jiro Yamada on 24 Apr 2019
ありがとうございます。
今のところだと、やはり右肩上がりのものが必要になるようです。

Sign in to comment.

Products


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!