光伝送シミュレーションを BASIC でやる。
なんだか非常に涙ぐましい努力に聞こえますが、
ちょっと試してみたい、というような場合にフリーソフトで
計算ができるのは、それなりに有用ではないかと思います。
ただ、もし、お仕事などでお使いになる場合は、プログラムの生産効率が
圧倒的に違いますので、エンジニアリング用のソフトをお使いに
なることを強くお勧めします。
ここでは、あくまで個人としての学習の一助となることを目的と
しています。
計算に使用している式などについては、
参考文献リスト
の光伝送シミュレーション関連の項にある文献を参照願います。
ここでは、まずその１として波長分散の影響のみを計算します。
シミュレーションは、いくつかのプログラムからなります。
これは、必要なものだけ個々に使用することが可能になること、
個々のプログラムの機能を明確にすること、
計算時間を分散することなどの理由によります。
全体の構成は以下のようになります。
それぞれのプログラムを十進 BASIC にてロードしては実行します。
各プログラムの計算結果は、プログラムファイルのあるディレクトリにある
data というフォルダの中に保存されます。
計算結果はテキストファイル形式ですので、
他のアプリケーションを用いてグラフ化するのも良いかと思います。
ただし、全データを生成すると11Mb以上になりますので、
ハードディスクの空き容量は十分確保してください。
なお、プログラムの実行については、
十進 BASIC のヘルプ、
「設定のヒント−＞アプリケーション実行向きの設定」にあるように、
プログラムファイルをダブルクリックにて実行できるようにすると便利です。
すべてのプログラムファイルはサンプルファイルに含まれています
（本ページのタイトル直下にダウンロードのリンクがあります）。
それぞれ、
-
prbs.bas
擬似ランダム符号系列（Pseudo-random bit sequence）を生成します。
任意の次数 n にて、
生成多項式、X^n + X^m + 1 = 0 （n, m はプログラム中で指定）の
ランダム符号列を n 個のシフトレジスタ（配列）を用いて生成します。
いわゆる、1,1,0,1,0,0,1,1,1,0,... というデータ生成です。
ここで、ちょっとした計算のコツですが、2^n-1 個のデータが生成されますので、
以降シミュレーションにおいて計算される FFT 演算の時間を短縮するため、
出来上がったファイルの最後にもうひとつのデータ、 0 を加えて、
（かつファイルの先頭にあるデータ数もひとつ増加させて）
2^n 個のデータとすると後の計算がずっと楽です。
プログラムに最初から追加しておけば良いという話もありますが．．．
-
NRZwaveform.bas
1,0,1,0,1,1,0,... という符号データ（PRBS符号列）より、振幅 2 (+/-1)、
オフセット 0 の NRZ 波形を生成します。
波形の立上がり（Tr）、立下り（Tf）部は正弦波を利用しています。
Tr、Tf 時間（0-100%）を一ビットピリオドの何％にするか設定できます。
また、クロスポイントも 10-90% の範囲で設定可能です。
波形の最初と最後の位相は連続していますので、
スペクトルを見る場合も FFT における窓関数は通常不要です。
-
MZmod.bas
マッハツェンダ干渉計光変調器のモデルです。
ひとつの電気変調波形の入力により、
最大振幅 1（絶対値）の変調された光電界波形を生成します。
マッハツェンダの２つの光路に対する変調効率を設定することにより、
チャープを設定することが出来ます。
生成されたデータは特に信号光波長に依存しないデータとなります。
-
dispersion.bas
FFTを利用して光スペクトラムを計算し、
それに、ファイバの位相に関する伝達関数β(ω)を
かけて波長分散の影響を分散スロープまで含めて計算するプログラムです。
入力は総分散値ですが、分散スロープとのバランスをとるため、
いったんプログラム内で距離に換算し、総分散を計算します。
またパワーは保存されますので、一切ロスのないモデルです。
正の分散および負の分散どちらでも
つかえますので分散補償としても使えます。
-
OE.bas
光電界波形から、パワーを計算し、
リファレンスフィルタ（ビットレートの 0.75 倍の帯域をもつ 4 次ベッセルフィルタ）
を通した波形と、
フィルタリング無し（理想 O/E 変換）の
波形を計算します。
実行時に、出力波形を平均 0.5 に規格化するか、
光ファイバ中のパワーに換算し、1[V/W] の規格化感度を使用するか
の選択ができます。
光ファイバ中のパワーへの換算は、
参考文献「Mathcadによる光システムの基礎」から引用しています。
-
eyediagram.bas
NRZwaveform.bas もしくは、OE.bas により生成された電気波形を、
ビットレートの整数倍で折り返してアイダイアグラムを作成します。
実際に計算を行う前に、
私から見たシミュレーションのポイントの説明を試みたいと思います。
シミュレーションは、光キャリア信号（一定パワーの光源）を
入力データ（1,1,1,1,0,0,0,1,1,0,...）にしたがってOn/Off強度変調した
信号に対して、波長分散の効果を計算します。
これは特に光に限ったものではなく、AM変調された信号の伝搬に関する
シミュレーションだともいえます。
さて、AM変調された信号を単純に数値計算しようとすると、
サンプリング定理を十分満足した以下のようなサンプル波形を処理する、
というのがまず最初に思いつくのではないでしょうか？
（私だけかもしれませんが...）
これを行うと、例えば、10Gb/s の信号の場合、
一ビットあたり 100ps になりますので、
1550nm（193.5THz）の光信号を変調したとすると、
サンプリング定理（最大周波数成分＝キャリア周波数の２倍でサンプル）
を満足しようとすれば
一ビットあたりのサンプル数は、
193.5THz × 2 × 100ps ＝ 38710サンプル／ビット
となります。
PRBS 2^7-1 という短いランダム符号を使ったとしても、
500万ポイントくらいのサンプルデータになってしまいます。
これは、パソコン＋BASICで計算するには若干きびしい見積もりです。
また、SPMなど時間軸上でおこる位相変調効果の計算は、
単純に瞬時強度をサンプルしたデータでは非常にややこしく（無理に？）なります。
ここで、導入されるのが複素数表示による変調信号の表現です。
具体的には、ωc をキャリア信号角周波数（光の角周波数）、
A(t) を上図の包絡線となる振幅をあらわすとし、
φを t=0 における初期位相角とすれば、
・・・・・(式1)
にてあらわされるAM変調信号 a(t) を、
・・・・・(式2)
と書き直し、a'(t) の実数部が実際の変調波形をあたえるとします。
なぜこのような置き換えが必要かというと、
「計算が非常に便利になるから」
といえます。
複素数を導入することは単にわかりにくくしているだけのように見えますが、
光伝送シミュレーションのような位相領域での演算が必要な計算に対しては、
劇的に計算が楽になります。
例えば、a(t)に位相変調 φm(t) を加えたいとします。式の上では、
・・・・・(式3)
・・・・・(式4)
となり、一見どちらも大差ないように思えます。
しかし、複素数表示（式４）の場合、
上の式を以下のように表すことができます。
・・・・・(式5)
つまり、元の信号（式２）に Exp(jφm(t)) をかけることで、
位相変調があらわせることを意味します。
これは、例えば元の信号が、a'(t) = f(t) のように具体的な形がわからなくても
位相に対する演算が可能であることを示します。
式１のような形の場合、
元の式が具体的にわからないと位相演算を行うことが出来ません。
この特性はサンプリングしたデータでもいえます。
たとえば、a(t)を時間 t0, t1, t2,... にてサンプルしたデータ、
a0, a1, a2,... に位相変調 φm(t) {t = t0, t1, t2,...}を加えようと
しても、いったんもとの式（式１や式３など）から計算しなおさなくては
なりません。
ところが、式２の形であらわされた a'(t)を同じように時間 t0, t1, t2,... にて
サンプルした複素データに対しては、
Exp(jφm(t)) を t0, t1, t2,...にてサンプルした値を個々にかけることにより
位相変調されたデータを得ることができます。
重ねて述べますが、各サンプルデータの実数部が実際の瞬時値になります。
また、式５を以下のように変形すると、
・・・・・(式6)
最後の Exp(jωt) は光のキャリア周波数成分であり、
第一項を、
・・・・・(式7)
とおくと、
・・・・・(式8)
このようにあらわすと、
フーリエ変換の周波数シフトの法則が適用できることがわかります。
つまり、a'(t) のフーリエ変換は、
f(t) のフーリエ変換を F(ω) とするなら、F(ω−ωc) にて与えられます。
逆に言うと、a'(t) のフーリエ変換を、
ωc＝0 として周波数ゼロの周りにシフトして計算を進めることが可能になります。
これを利用し、
ここでのシミュレーションでは、FFT のサンプルを ωc＝0 とし、式７であらわされる
位相情報を含む包絡線成分をサンプルすることにより計算を進めます。
サンプルされたデータを、
ここでは光サンプルデータと呼ぶことにします。
これにより、包絡線の変化は一ビット内でゆっくりと変化しますので（前提です）、
サンプル数も、一ビットあたり 64 としています。
先にあげた38710サンプル／ビットに比べると 600 分の一の情報量ですみます。
分散などの計算もすべてキャリア周波数をゼロとして、
キャリア周波数からの偏差 Δω にて計算されることになります。
この光サンプルデータから光の波形を得るには、
各サンプルデータの大きさをとり包絡線を作成すればよいことにご注意ください。
実数部をとるのはあくまでキャリア周波数成分まで含んだ表現の場合です。
複素数による表現にて一点注意が必要なのは、
信号同士の掛け算は注意が必要ということです。
例として、初期位相のことなる以下の二つの信号を例に取ります。
・・・・・(式9)
この二つの信号の積は、
・・・・・(式10)
が正しい解ですが、式９を複素数表示を用いて、
・・・・・(式11)
とあらわして、x'(t)×y'(t) を計算しても式１０にならないことは明らかです。
このことは、参考文献にあげている、Yariv の光エレクトロニクスの
冒頭にて述べられています。
変なたとえですが、パワーをあらわす dBm という単位は、ロスや利得などの
計算が楽になるように導入されましたが、パワーの和は dBm でやっては
いけないというのと似ていると思います（1dBm＋2dBm は 3dBm とはなりません）。
サンプルプログラムを例にとり、シミュレーションの実行をおってみます。
主なシミュレーションの条件は、
- PRBS 2^9-1
- ゼロチャープ
- 総分散量は実行時に入力
- 分散スロープの影響まで含む
です。
まず、prbs.bas によりランダムデータを生成します。
結果は、data フォルダに prbs.txt として保存されます。
後で行う FFT の計算を楽に（早く）するため、
この prbs.txt ファイルをテキストエディタで開いて、
データの最後に 0 を加え総データ数を 2^n 個にします。
具体的には、ファイルの先頭にある、
1
REAL
511
この３行目が総データ数ですので、これを１加えて 512 とし、
ファイルの最後に
0
を一行加えます。
サンプルファイルでは、この結果を prbs_p0.txt として保存してあります。
次に、NRZwaveform.bas にて変調用の NRZ 波形を生成します。
このとき、入力として先ほど述べた prbs_p0.txt を使用しています。
結果は data フォルダに NRZwaveform.txt として保存されます。
ご参考までに、
この結果をエクセルを使用して最初の 100bit 程度の波形をグラフ化したものが
以下となります。
eyediagram.bas により NRZwaveform を表示したものが以下です。
この例では、特にＸ軸・Ｙ軸のラベルの表示を行っていません。
必要でしたら適当にプログラムを変更するなどして加えてください。
この生成された変調用波形を用いて、
MZmod.bas により光信号を得ます。
正確には、前のシミュレーションの説明の章にて述べた、
位相情報を含む包絡線サンプルデータを得ます。
ここで使用している変調器のモデルは、
以下のように入力された光パワーが各導波路（１，２）に
等しいパワーで導波されるとして各導波路に振幅 1/2 を与え、
さらに各導波路長は仮想的にゼロとし、
変調電圧に比例した位相シフト−φ１および＋φ２を受けるとします。
また、分岐・合成による位相変化も無しとしています。
位相シフトはシミュレーションの説明の章で述べたように、
Exp(jφ) を信号にかけることで行なえます。
φ１＝φ２（チャープフリー）の場合の、
出力の光電界強度 Ｅ および光パワー Ｐ と総位相差 φ１＋φ２ の
関係を示します。Ｅ および Ｐ は位相差に対して余弦波関数となります。
今回の場合位相差は変調電圧に比例しますので、結局、変調電圧に対して
余弦波（正弦波）型の出力パワー関数が得られます。
（正確には Ｅ は透過光強度の振幅を表すべきで、絶対値を用いてあらわされるべきですが、
以下の図では余弦波であることが理解しやすいよう位相情報をふくんだ＋/−の値で図示しています。）
なお、位相差 π をあたえる変調電圧を Vπ といいます。
変調電圧に対する各位相シフト量の効率（結合係数）は、プログラム中で設定します。
この効率と変調電圧ゼロにおける初期位相差（バイアス点）設定により
チャープを設定します。計算結果は、mzModOut.txt として保存されます。
今回の場合、とりあえずチャープフリーになるように設定してありますので、
位相が固定された複素数データ（虚数部＝ゼロ）が生成されます。
チャープ量の確認については、各データ間の
位相差をサンプル時間間隔で割ることにより
キャリア周波数からの周波数シフト量が得られます。
さて、生成されたデータ mzModOut を、
OE.bas にて電気信号（実数データ）に変換し、
平均を 0.5 に規格化した結果（idealOE.txt）をアイダイアグラムに
したものが以下です。
次に、波長分散計算をおこないます。
まず、光サンプル波形の複素 FFT により得られた n 個の係数を、
X1〜Xn (n：データ数）とし、
DFT の周期性を利用して X(n/2+1)より上のデータを
マイナス周波数側にシフト（移動）します。
DFT の周期性については、
本サイトにある「離散フーリエ変換（DFT）と高速フーリエ変換（FFT）」
の例が参考になるかと思います。
得られたスペクトラムに、各周波数成分に対応するファイバ伝送による
位相シフト量を加え、分散の影響を計算します。
以下に計算の様子を示します。
ここで、波長分散の影響を計算するためには、
位相シフト量は定数ではなく角周波数 ω に依存した値として扱う必要があります。
しかし、そこに踏み込む前に、まず位相角の方向を以下のように確認します。
入力点（原点）から距離 L[m] だけ離れた点での位相は、
位相定数をβ[rad/m]とすれば、原点での位相にくらべて βL[rad] だけ前の位相となります。
したがって、シミュレーションの説明のところでも若干ふれたように、
ある時間 t における入力点の信号に、
exp(-βL) をかけることで距離 L だけ離れた点での時間 t における位相を知ることができます。
そして今、βはωに依存した関数であることがわかっていますが、とりあえず定数だと仮定します。
すると、exp(-βL)も定数となりますので、
これは周波数領域でもそのまま exp(-βL) をかけることになり、
このとき位相角φは -βL[rad] となります。
さて、ここで、あらためて位相定数βを周波数領域にて各周波数ωの関数であるとし、
βをキャリア角周波数 ωc の周りにテーラー展開した式を以下のようにおきます。
・・・・・(式13)
β0、β1、．．．はそれぞれωcにおけるテーラー級数の
係数ですが、以下のような意味を持ちます。
- β0：周波数に依存しない位相オフセット値
- β1：群遅延時間
- β2：分散
- β3：高次の分散
これは、群遅延が τg＝−dφ/dω＝d(βL)/dω となり、
さらに群遅延を ω にて微分したものが分散になることから理解できます。
分散の影響を計算する場合は、波形ひずみが興味の対象ですので、
位相オフセット値および群遅延時間に関しては無視します
（波形ひずみに影響を与えないため）。
今、光スペクトラムは ωc がゼロにシフトされていますので、
ω−ωc＝Δω としてスペクトラムの各周波数成分に対応する
位相定数を計算します。
以下に、波長 1550[nm] において分散 16[ps/(nm km)] および
分散スロープ 0.06[ps/(nm^2 km)] としたときの、
式１３を利用した波長−群遅延特性を示します。
縦軸は群遅延 [ns/km]、横軸は波長 [nm] です。
1550[nm] を基準としていますので群遅延がその点でゼロとなっています。
計算に使用したプログラムは checkDispersion.bas です。
単位変換の例として参照いただけるのではないかと思います。
同じプログラムにて計算された、波長−分散特性の計算例です。
縦軸は分散 [ps/(nm km)]、横軸は波長 [nm] です。
さて、実際の波長分散計算をつづけます。
マッハツェンダ変調器の出力波形 mzModOut.txt を
入力として、dispersion.bas を総分散値 1600[ps/nm] とし実行します。
得られた結果は光波形ですので、
OE.bas を用いて平均 0.5 に規格化された電気信号に変換し、
アイダイアグラムを作成した結果が以下です（idealOE.txt）。
リファレンスフィルタを通した波形は以下のようになります。
