HB9CV向けアンテナシミュレータ

１． はじめに

八木宇田アンテナのシミュレータはいろいろありますので、いまさら新しく作る必要は無いのかもしれません。しかし、八木宇田アンテナの動作の理解を深めるためと、自分独自のプログラムがあればいろいろと改良ができて便利であろうという理由でアンテナシミュレータを作ってみることにしました。
ただし、数値計算から作っていると時間がかかり、アンテナ以外のことをいろいろやらなければなりません。そこで、信号処理などの計算に用いられるMATLABというツールを使ってシミュレータを作ることにしました。

２． MATLAB

MATLABはMathworks社のシミュレーション用プログラミング環境です。MATLAB独自のプログラミング言語が組み込まれており、M-スクリプトとかM-言語とか呼ばれています。MATLABはインタープリタ言語とその動作環境を含んだようなプログラムで、昔のパソコンに付属していたBASICインタープリタのような使い方をします。ただ、数値計算、特に行列計算に関する機能は非常に強力で、信号処理や制御、金融工学などの分野で標準的に用いられています。
言語処理系として見た時の特徴はAPLとC言語の便利な部分を抜き出して合わせたような仕様の言語であるということです。この両者の言語の経験者は何の苦も無くMATLABでプログラミングできるでしょう。ただし、APLなどというプログラミング言語はあまり使われていませんので、大部分の人はAPLのような動きをする部分で苦労するのではないかと思います。C言語のような手続き型言語ユーザーが最も苦労するのは、いかにしてループ処理を書かずにプログラミングするかに尽きるでしょう。また、行列の数学的知識も必要です(しかし、これが無いと信号処理の計算自体が理解不のですね)。
このように人によってはとっつきにくい言語ですが、かなり高水準なプログラミング言語ですので、非常に見通しの良いプログラムを書くことができます。

３． アンテナシミュレータの概要

今回作成したアンテナシミュレータは起電力法という方法で計算しています。各エレメントに正弦波形状の電流が流れるものと仮定して計算します。起電力法については文献「八木アンテナを作ろう―電脳設計ソフト　ＹＳＩＭで作る八木アンテナ (」ISBN:4789810690)」等を参照してください。
ここでは簡単にですが、計算手順を書きます。次節にソースコードを挙げますので、詳しくはそちらを参照してください。
(1)起電力法
今回用いた方法は起電力法という方法で、各エレメントに正弦波電流が流れているものと仮定して計算します。エレメント長が波長の1/2に近いとき、この仮定はそれなりに良い近似であるようです。
(2)エレメント自己インピーダンスの計算

エレメントの自己インピーダンスを求めます。ここではこのインピーダンスはそのエレメントが単独に存在した場合に、導体の中心から給電した時のインピーダンスをこう呼ぶことにします。

よく、ダイポールアンテナの自由空間での入力インピーダンスは約73Ωだと聞きますが、これはここの計算により求めることができます。
自己インピーダンスの計算は、まず円筒導体に加わる電力を求めることから考えます。起電力法を用いるため、導体上の電流は正弦波状に分布すると仮定します(図1)

。

図1．導体に流れる電流

したがって、導体の長軸方向をz軸とすると、導体に流れる電流は式(1)のように計算できます。
                                            (1)
ここで、Ioは電流の最大値、?は位相定数、zは導体端よりの距離です。位相定数は距離をかけると位相となる数値で、波長を?とすると式(2)のようになります。
                                                       (2)
この位相定数は距離から位相情報を算出する時に使用し、後述のビームパターン計算にも用います。
次に、円柱導体表面の起電力を求めます。電磁気学より位置zでの微小区間における導体表面の電界は次の式(3)のようになります。
                          (3)
ここで、Ex(z)は導体の中心から距離zの位置の電界、rは導体の半径、r1、r2は導体端からp点までの距離、r3は導体中心からp点までの距離です(図2)。

図2.起電力を求めるためのパラメータ

そこで、式(4)を円筒導体全体に渡って積分すると、導体中心から見た円筒導体表面の起電力を求めることができます(式4)。
                  (4)
これで起電力と電流の式が決定できましたので、導体に加わる電力が分かります。電流が正弦波状に分布するので、導体の中心が電流最大値となるはずです。また電位中心での電位をEとすると、電力Pは式(5)となります。
                            (5)
円筒導体の自己インピーダンスは式(6)で定義できます。
                                      (6)
式(6)に式(2)、式(5)を代入すると式(7)になります。
                              (7)
式(7)の分子と分母にIoがあり、Ioは0でないため、式(7)中の分子分母をIoの2乗で割ります。
                               (8)
式(4)の積分に式(3)中のIoは定数としてかけられているだけであるため、式(8)に式(4)を代入するとIoが完全に消去され、式(9)となります。
                           (9)
この式(9)を計算することにより、自己インピーダンスを求めることができます。
ただし、r1,r2,r3,はzの関数であるため、
                                                                           (10)
として計算します。
プログラムでfPz.mが式(9)の分子と式(10)の計算に相当します。MATLAB上での計算はリストリスト(dipoleimp.m)のようになります。この計算結果はZ = 73.1280 +41.9163i[?]となり、73Ω近辺の値となります。
MATLAB関数quad8は数値積分を求る関数です。q = quad('fun',a,b,tol,trace,P1,P2,...) は、指定した関数 G = fun(X,P1,P2,...) に係数 P1, P2, ... を直接渡します。tol や trace にデフォルト値を使っているため、空行列を渡しています。
(3)エレメント相互インピーダンスの計算
相互インピーダンスは自己インピーダンスを求めるときに導体表面の起電力を求めたところを、相手側のエレメントの起電力に変更することで計算できます。
詳細は工事中！！
自己インピーダンスの計算との変更点は僅かです。リストfmuPz.mが電位/Ioの部分を計算する関数で、quad8で積分しています。
(4)エレメント電流の計算
エレメントの自己インピーダンス、相互インピーダンスが求められると、電気回路に置き換えた回路解析が可能です。
詳細は工事中！！
(5)ビームパターンの計算
詳細は工事中！！
追加で書き足してみていますが、まだまだ未完成です。未完成なのにWeb上なのは、これ以上遅れると忘れてしまいそうだからです。なので、まだまだ不備が多いです。ただ、だいたいこんな感じになるということですので、こうした方が分かりやすいとかであれば、ご連絡ください。
(5-1)平面波の到達
電波をどちらの方向に放射するのかを求めるには、どの地点で見た場合かを決めなければなりません。アンテナのすぐ近くで見た場合とずっと離れた場合では違った広がりになります。無線通信に用いる場合は近い場合よりもずっと離れている場合が多く、無限遠点で計算すると計算が簡単であるため、無限遠点を基準にビームパターンを計算することとします。無限遠点を基準に考えると、平面波と考えることができます。平面波とは波面が平面である波です。言い換えると、その平面内には同一時刻に波面が到達します。

図３． 平面波

立体的に描くと分かり難いので、次元を減らして2次元の絵で描いてみることにします。

図４． 次元を減らして描いた平面波

こうするとかなり分かり易くなります。図中の線(実は平面)の上に同時に波が到達するというのが平面波です。
シミュレーションをする場合、エレメントを分割して計算する必要があります。そこで平面波が分割したエレメントにどのように到達するのか考えます。

図５． 分割したエレメントに平面波が到達

平面波がエレメントと直角方向からやって来た場合、分割したエレメントには同時に電波が到達します。絵で見ているとエレメントの位置によって到達時間が違うように思えますが、そうならないのが平面波なのです。電球の光は広がっていて、その光でできる影は大きくなりますが、太陽だとそうなりません。太陽は地球から遠いため、その光はほぼ平面波だからです。無線の場合も相手が十分に遠いので、平面波と近似しても問題ありません。
それでは、次は斜めから電波が到達した場合です。

図6.入射角θの平面波の到達

2つの離れた地点に斜めから平面波が到達した場合、当然その2地点に平面波の到達する時間が変わります。今図のようにx離れた地点に入射角θ方向から平面波が到達したとします。すると波の伝わる経路はlだけ長くなります。この長さlは次の式のように計算できます。
したがって、この経路差だけ波が遅れて到達します。一方、電波は繰り返して到達するので、時間が遅れるというのは位相が遅れたのと同じ事になります。また、位相の遅れとして表現した方が簡単になります。そこで、単位距離の経路差によってどれだけ位相がずれるかを表す位相定数βを使うことにします。
となります。つまり1波長の経路さは2?の位相変化に相当するということです。
したがって、式より入射角θの場合、x離れた位置での位相遅れは次の式のようになります。
これに式を代入すると次の式のようになります。
一方、電波は周期的に繰り返される波です。このため、無変調波の電界は次の式のように複素関数で表現することができます。ここでjは虚数単位とします。
これに位相遅れが加わると
となります。ここではx点の電界、は振幅を示す定数です。あまり文字が多いと分かりにくくなるため、以後とすることにします。
(5-2)特定方向への放射
前節で分割したエレメントに平面波が到達する場合の位相遅れ量の計算ができるようになりました。送信と受信はほぼ逆の動きをすると考えられますので(受信した場合の電流分布は異なるが、平面波しか扱わなければ同じ)、前節の方法で送信を考えることができます。
話を簡単にするため、2つの点波源から?方向に電波を放射する方法を考えてみます。

図７．2つの点波源からの放射

?方向に強く電波を放射するためのx離れた点波源の位相差は受信の時と同様に考えると
となります。このため、次の式のように点波源を励振すると?方向に波が出ます。
複素関数で書くと便利なところは、位相の足し算が掛け算で表現できるところです。上の式は次の式のように書くことができます。
ところで、オイラーの定理より
となるので、式は
と書くことができます。つまり、普通の複素数の掛け算で位相のずれを計算できることになります。また、式と式は複素共役の関係にありますので、式の共役複素数を求める方法によっても計算することができます。
(5-3)複数の等方性アンテナからの放射
計算方法の理解のため、前節までの内容をプログラムで組んでみます。
ex1.m
このプログラムは5個の当方性アンテナを使い30度方向にビームを振った例です。30度方向に振ると、180度−30度＝150度にもローブが出てしまいました。

図8．プログラム実行結果

(6)ゲインの計算
詳細は工事中！！

４． プログラム本体

そのうちにプログラムのドキュメントを書く予定です。ソースリストはアップしますので、読んで見てください。また、まだ正確かどうか自信がありません。誤り等をご指摘いただければ幸いです。
また、MATLABは非常に高価な処理系ですので、誰でも動作させるというわけには行きません。そこで、アマチュア向けにはGNU Octave(http://www.octave.org/)などの処理系が良いのでしょうが、まだOctaveでの確認はできていません。
プログラム本体(ソースリスト)
エレメント自己インピーダンス計算用関数
エレメント相互インピーダンス計算用関数
ビームパターン表示ルーチン
エレメント電流表示ルーチン
テストドライバプログラム

５．実行結果

図1 .50MHz 8エレエレメント電流

図2.50MHz8エレビームパターン

図3. 50MHz3エレHB9CVエレメント電流

図4. 50MHz3エレHB9CVビームパターン

図5. 14MHz3エレエレメント電流

図6. 14MHz3エレビームパターン