MATLABのTips 伝達関数行列の効率的なコーディング
伝達関数行列とは
今回は伝達関数行列のコーディングについて紹介したいと思います。
ちなみに伝達関数行列は、各周波数における伝達関数を入力×応答数分に配置した行列です。
音場制御や伝達経路解析（TPA）などで良く使用します。
ここで、伝達関数行列Gは入力ベクトルqと応答ベクトルpとの間に下式の関係があります。
伝達関数行列の要素は下のようになります。
※要素内も角周波数ωの関数となりますが、見にくくなるので省略しています。
ここでもし学生さんに「G11(ω)などは何なのでしょうか？」と問われたら、「伝達関数です。」という回答になってしまいます。
この辺で頭が混乱しだすと思います。
というのも、伝達関数自体は測定器などで簡単に取得できますが、伝達関数行列の構成方法が良く理解できないからです。
ここについては”式みたらわかるじゃん”的な発想で省略されているので、独学で学習している方には優しくない点だと思います。
ここでは超優しく解説したいと思います。
実験や解析で取得した伝達関数（FRF）は下記のようなイメージかと思います。これが行列の1要素となります。
つまり、上記FRFがG1,1を指す場合、伝達関数行列の1行1列目に配列され、且つ角周波数ωに応じて値を変える必要があります。このときFRFは「実数＋i×虚数」で表現する必要があります。
測定機によっては、FRFは振幅と位相に分かれたり、デシベルなどの単位でエクスポートされることがありますが、伝達行列にする際には単位を変換して「実数＋i×虚数」にする必要があります。
例えば、音場制御で制御用のスピーカーが10個でマイクロフォンを100本で測定したとき、Inputは10でOutputが100なので、100×10の伝達関数行列になります。この伝達関数行列は周波数に応じて中身の伝達関数を変更する必要があります。
3階のテンソル or セル
セル
だいたいの学生が伝達関数行列のような行列を構成する場合、ぐちゃぐちゃなコードを書くか、cell関数を使うのではないかと推察します。
cell関数も悪い手ではないです。（私も学生のときはcell関数で構成していました。）
ちなみに、cell関数では、cell(N,M)と入力すると、N行M列のcellが作成され、その中にデータを格納できます。
cellは箱のようなものなので、任意のデータ形式を格納できます。
下記では、1行1列にベクトルを格納し、1行2列にストラクチャー構造のデータを格納しました。
> >ce=cell(2,3)
ce =
[] [] []
[] [] []
ce{1,1}=[1;2;3];
> > a.c=1;
> > ce{1,2}=a
ce =
[3x1 double] [1x1 struct] []
[] [] []
上記のようにcellの中にFRFを格納してωに対応した値を呼び出すというのでも、問題ありません。
ただし、格納したデータを呼び出すのが面倒で、呼び出すためのコードを書く必要があります。
そのため、私はあまり推奨しません。
具体的には、下記例のようにbというcellを作成し、各cell内にaというベクトルを格納します。
そのときに、全cellのa(5)を抽出したい場合を考えてみてください。
この場合for分を使うしかありません。（自分でトライ＆エラーしてみてください。）
clear all
a=1:10;
b=cell(2,3);
for ii1=1:size(b,1)
for ii2=1:size(b,2)
b{ii1,ii2}=a;
end
end
% b{1:2,1:3}(1)
3階のテンソル
「テンソルって？」と質問された場合、詳しくないので答えにくいのですが、下図が参考になると思います。
この図をみると、そもそも伝達関数行列って3階のテンソルじゃない？と気が付いた方がいると思います。
そうです。伝達関数行列は3階のテンソルです。
なので、最初から3階のテンソルとして、データ構成を考えた方が楽です。
「じゃあなんで最初から3階のテンソルとして扱わないの？」と疑問に感じると思います。
理由は、zerosやonesには、zeros(M,N)のように行列を扱えることの説明はありますが、3階のテンソルを扱える説明が弱いからです。
下記がzerosやones関数の説明ですが、「ZEROS(M,N,P,...) または .....」というのが説明になります。学問自体でテンソル計算をする場合は、学生でも気が付くのだと思いますが、振動音響工学は行列計算がメインなので、3階のテンソルで表現できる！ということになかなか気が付きません。
一般的にコードを書く場合は、zeros関数などで、データのメモリをあらかじめ確保して、zerosで作成した行列にデータを格納します。
つまり、zerosで3階のテンソルが扱えるということに気が付けなければ、"3階のテンソルを使う"という考えすら思いつかないのです。（思いつかない人が多いのです。）
>> help zeros
ZEROS 零の配列
ZEROS(N) は、N×N の 0 の行列です。
ZEROS(M,N) または ZEROS([M,N]) は、M×N の 0 の行列です。
ZEROS(M,N,P,...) または ZEROS([M N P ...]) は、M×N×P×... の 0 の
配列です。
ZEROS(SIZE(A)) は、A と同じサイズで、すべて 0 です。
引数なしの ZEROS は、スカラの 0 になります。
ZEROS(M,N,...,CLASSNAME) または ZEROS([M,N,...],CLASSNAME) は、
クラス CLASSNAME の 0 の M×N×... の配列です。
注意 : サイズの入力 M,N,P... は、非負の整数である必要があります。
非負の整数は、0 として扱われます。
例:
x = zeros(2,3,'int8');
参考 eye, ones.
Overloaded methods:
darray/zeros
Reference page in Help browser
doc zeros
>> help ones
ONES すべての要素が 1 の配列の作成
ONES(N) は、N×N の要素がすべて 1 の行列を出力します。
ONES(M,N) または ONES([M,N]) は、要素がすべて 1 の M×N の行列を
出力します。
ONES(M,N,P,...) または ONES([M N P ...]) は、要素がすべて 1 の
M×N×P×... の配列です。
ONES(SIZE(A)) は、A と同じサイズですべての要素が 1 の配列です。
引数なしの ONES は、スカラ 1 です。
ONES(M,N,...,CLASSNAME) または ONES([M,N,...],CLASSNAME) は、クラス
CLASSNAME の 1 からなる M×N×... の配列です。
注意 : サイズの入力 M,N,P... は、非負の整数である必要があります。
非負の整数は、0 として扱われます。
例:
x = ones(2,3,'int8');
参考 eye, zeros.
Overloaded methods:
darray/ones
Reference page in Help browser
doc ones
例えば、zeros(2,5,3)のように入力すると3階のテンソルが作成できます。
>> zeros(2,5,3)
ans(:,:,1) =
0 0 0 0 0
0 0 0 0 0
ans(:,:,2) =
0 0 0 0 0
0 0 0 0 0
ans(:,:,3) =
0 0 0 0 0
0 0 0 0 0
3階のテンソルとして、Inputは2でOutputが3の伝達関数行列を簡易的に作成してみましょう。
各FRFは適当なベクトルを用意します。なお、角周波数ωとFRFベクトルの長さは一致するので、ω(n)のときの伝達関数の値はFRF1_1(n)のようになります。したがって、ω(n)のときの伝達関数行列はG(:,:,n)となります。
ここで、協調したいことは、構成した伝達関数行列を行列として取り出すのが用意だということです。
G(:,:,n)と一行のコードで取り出せるということが複雑なコードを書く上で重要となります。
clear all
omega=0:9;
G=zeros(3,2,10);
FRF1_1=[1 2 3 4 5 6 7 8 9 10]*1;
FRF2_1=[1 2 3 4 5 6 7 8 9 10]*2;
FRF3_1=[1 2 3 4 5 6 7 8 9 10]*3;
FRF1_2=[1 2 3 4 5 6 7 8 9 10]*1*5;
FRF2_2=[1 2 3 4 5 6 7 8 9 10]*2*5;
FRF3_2=[1 2 3 4 5 6 7 8 9 10]*3*5;
G(1,1,:)=FRF1_1;G(1,2,:)=FRF1_2;
G(2,1,:)=FRF2_1;G(2,2,:)=FRF2_2;
G(3,1,:)=FRF3_1;G(3,2,:)=FRF3_2;
% % ω=0
G(:,:,1)
今回はこのへんでGood luck
ここでは扱っていない機械エンジニアネタもありますので、下記サイトもよろしくお願いいたします。
Discussion