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

以下のコードで得られるボード線図の伝達関数を求めたいです。

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)のボード線図となる伝達関数を求めることができるのでしょうか?

Answer 1

[Hiroumi Mita]

掲題の問題について

以下のようにシステム同定の手法でモデルを推定しますと

一応、ほぼにたようなボード線図を書くことはできます。

しかし、得られた伝達関数は複素係数のものになり

現実性は乏しく、無理やり数値計算で合わせこんだ印象があります。

もともと推定したいモデルのボード線図は

周波数帯域において、ゲイン特性が一定数(0dB)、位相は

周波数が上がるとじりじり上がるもので、現実的に

このようなモデルがあるのかどうか、必要性はあるのか、

数式表現できるものなのか大いに疑問があります。

この問題の前提条件を再確認されることをお勧めします。

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})%ボード線図

[Jiro Yamada]

お返事が遅くなって申し訳ありません。

回答ありがとうございます。

恥ずかしながら伝達関数について理解が足りなく、位相が上がっていく特性がおかしいことは指摘されて初めて知りました。

前提条件を見直してみます。

[Jiro Yamada]

度々すみません。

位相が下がっていくこちらならどうでしょうか。

同じコードで試したのですが、なかなかにかけ離れたものが得られてしまいました。

どこを調節すればよいかなども教えていただけると幸いです。

たくさんお願いして申し訳ありませんが、よろしくお願いします。

%----------パラメータ設定-----------%
 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})%ボード線図 

Answer 2

[Hiroumi Mita]

ボード線図において、例えば次のような特徴を持つものは

無駄時間が該当します。

#1. ゲイン特性が周波数帯域において一様0[dB]

#2. 位相特性は周波数が高くなると遅れる

このサンプルが無駄時間1[s]のボード線図です。

s=tf('s')

sys=exp(-s)

figure;bode(sys)

無駄時間は入力が入っても一定時間（無駄時間）無反応なものを

言います。

伝達関数は

sys=exp(-T*s)

T:無駄時間

です。

もし、このシステムを対象にするなら、出力が何秒遅れて入力に反応したかを

調べればよく、ボード線図から推定するのは意味がありません。

何らかの動特性を持つシステムは、無駄時間のみ表現されるものはないと思われます。

動特性（時間による変化）＋無駄時間が一般的です。

基本的に、ボード線図からモデルを推定するとき、周波数の変化に対する

ゲイン特性の変化、位相の変化を必要とします。

そこで、示されたボード線図は、何らかの実験の不手際により、

ゲイン特性が正しく取れなかったか？、元の仮定では動特性を無視したのか？

そんなことが考えられるわけです。

動特性は示された系に無いのか？もし、本当にないのならそれはボード線図で扱う領域ではないかもしれません。

そのあたりを再考いただければと思います。

[Jiro Yamada]

解説までしていただきありがとうございます。

Simulinkで変位センサの信号を通すデジタルフィルタを作りたくて伝達関数の推定を行っております。

希望のボード線図は実験で得られたものではなく、数値計算を基に求めたものになります(最初のコードの13行目にあるF1=cから)。

教えていただいたことを参考に考えてみます。

Answer 3

[Shoumei]

位相が右上がりではなく右下がりで良ければこんな特性のフィルタはありますが。

ギターのエフェクトでよく使われるフェイザーのオールパスフィルタです。

[Jiro Yamada]

ありがとうございます。

今のところだと、やはり右肩上がりのものが必要になるようです。