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

23 次查看（过去 30 天）
Jiro Yamada
Jiro Yamada 2019-4-5
评论： Jiro Yamada 2019-4-24
以下のコードで得られるボード線図の伝達関数を求めたいです。
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)のボード線図となる伝達関数を求めることができるのでしょうか?

请先登录，再回答此问题。

采纳的回答

Hiroumi Mita
Hiroumi Mita 2019-4-8
掲題の問題について
以下のようにシステム同定の手法でモデルを推定しますと
一応、ほぼにたようなボード線図を書くことはできます。
しかし、得られた伝達関数は複素係数のものになり
現実性は乏しく、無理やり数値計算で合わせこんだ印象があります。
もともと推定したいモデルのボード線図は
周波数帯域において、ゲイン特性が一定数(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 个评论
Jiro Yamada
Jiro Yamada 2019-4-17
お返事が遅くなって申し訳ありません。
回答ありがとうございます。
恥ずかしながら伝達関数について理解が足りなく、位相が上がっていく特性がおかしいことは指摘されて初めて知りました。
前提条件を見直してみます。
Jiro Yamada
Jiro Yamada 2019-4-17
度々すみません。
位相が下がっていくこちらならどうでしょうか。
同じコードで試したのですが、なかなかにかけ離れたものが得られてしまいました。
どこを調節すればよいかなども教えていただけると幸いです。
たくさんお願いして申し訳ありませんが、よろしくお願いします。
%----------パラメータ設定-----------%
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})%ボード線図

请先登录，再进行评论。

更多回答（2 个）

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

请先登录，再进行评论。

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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 データの解析 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!