【Serdes】Dual Dirac MethodをPythonで実装

Dual Dirac MethodをPythonで実装してみた。

Dual Dirac Methodの技術参考資料 (Keysight)
https://www.keysight.com/upload/cmc_upload/All/dualdirac1.pdf

1. Jitter Dataの準備

解析するモチーフをまず定義する。
Data Rate: 10Gbps (1UI=100psec)
Rj: 1ps RMS
Dj: 14ps PP (±7ps)のSj
Sample Size: 1MUI

Jitter Dataのデータは、各行に各UI毎のJitter値として、text file dataを作成。
RjとDjそれぞれのデータを作成してから、それらの和を取ることでTjのデータを作成することにした。

1.1 Rjの生成プログラム

乱数でRjデータを発生させ、rj.csvに保存。

import numpy as np
import csv

#Sample size
intSize = 1000000

#Jitter RMS = 1ps
fltJitterRms = 1e-12

#Generate random pattern
arrData = np.random.normal(0, fltJitterRms, intSize)

#Output Data
arrWriteData = []
for i,d in enumerate(arrData):
        arrWriteData.append([])
        arrWriteData[i].append(d)

with open('rj.csv', 'w') as f:
        wr = csv.writer(f)
        wr.writerows(arrWriteData)

1.2 Dj(Sj)の生成プログラム

14ps-pp, 周波数100MHzの正弦波でSjを発生し、sj.csvに保存。

import numpy as np
import csv

#Sample size 
intSize = 1000000
#Time step : 1UI = 100ps
fltStepTime = 100e-12

#Noise frequency = 101MHz
fltNoiseFreq = 101e6
#Noise Amp = 14ps-pp
fltNoiseAmp = 14e-12 / 2

#Generate sin pattern
arrTime = np.arange(0, intSize*fltStepTime, fltStepTime)
arrData = []
for t in arrTime:
	d = fltNoiseAmp*np.sin(2*np.pi*fltNoiseFreq*t)
	arrData.append(d)

#Output Data
arrWriteData = []
for i,d in enumerate(arrData):
	arrWriteData.append([])
	arrWriteData[i].append(d)

with open('sj.csv', 'w') as f:
	wr = csv.writer(f)
	wr.writerows(arrWriteData)

1.3 RjとDjの合成

rj.csv, sj.csvをコマンド因数で受け取り、2つのジッタの和を取り、rjsj.csvに保存。

import numpy as np
import csv 
import sys

args = sys.argv
path = []
path.append(args[1])
path.append(args[2])

arrData = []

#Read path[0] data
with open(path[0], 'r') as f:
	rr = csv.reader(f)
	for rd in rr:
		arrData.append(float(rd[0]))

#Add path[1] data
with open(path[1], 'r') as f:
	rr = csv.reader(f)
	for i,rd in enumerate(rr):
		arrData[i] += float(rd[0])

#Output data
arrWriteData = []
for i,d in enumerate(arrData):
	arrWriteData.append([])
	arrWriteData[i].append(d)

with open('rjsj.csv', 'w') as f:
	wr = csv.writer(f)
	wr.writerows(arrWriteData)

1.4 生成したジッタデータの確認

結果確認用にヒストグラムプロットのプログラムも作った。

import numpy as np
import csv
import sys
import matplotlib.pyplot as plt

args = sys.argv
path = args[1]
arrData = []

#Read data file
with open(path, 'r') as f:
        rr = csv.reader(f)
        for rd in rr:
                arrData.append(float(rd[0]))

#Plot Histgram
plt.figure()
plt.hist(arrData, bins=100, normed=True)
plt.show()

1.1で作ったrj.csvヒストグラム。ちゃんと正規分布が得られている。

f:id:kent_s:20180817140053p:plain
rj.csv

1.2で作ったsj.csvヒストグラム。正弦波による特徴的なジッタ分布が得られてる。

f:id:kent_s:20180817140241p:plain
sj.csv

1.3で作ったrjsj.csvの分布。正弦波に正規分布が重畳されているのが分かる。

f:id:kent_s:20180817140407p:plain
rjsj.csv

2. Dual Dirac Methodの実装

2.1 BERプロットの作成

1.3で作ったrjsj.csvを使ってBERプロットを作成。
1UIを1000等分にスライスし、ジッタがx番目の各スライスに含まれる確率をJ(x)と定義。
次にJ(x)の積分により、BER(x)を求めています。
出力はber.csvへのファイル出力です。

import numpy as np
import csv 
import sys
import matplotlib.pyplot as plt

fltUi = 100e-12
intRange = 1000

args = sys.argv
path = args[1]
arrJx = np.zeros(intRange, dtype='float')
intSize = sum(1 for line in open(path))

# Read Jitter Data
with open(path, 'r') as f:
	rr = csv.reader(f)
	for rd in rr:
		d = float(rd[0])/fltUi*intRange
		idx = int(round(d))
		arrJx[idx] += 1/float(intSize)

# Cal BER
arrBer = np.zeros(intRange, dtype='float')

for i in range(intRange/2, 0, -1):
	arrBer[i] = arrBer[i+1] + arrJx[i]

for i in range(-intRange/2, 0, 1):
	arrBer[i] = arrBer[i-1] + arrJx[i]

arrBer[0] = (arrBer[-1]+arrBer[1])/2 + arrJx[0]

#Generate Time axis
arrT = np.arange(0, fltUi, fltUi/intRange)

#Plot BER
plt.figure()
plt.plot(arrT, arrBer)
plt.yscale('log')
plt.show()

#Write BER
arrWriteData = []
for i,d in enumerate(arrBer):
	arrWriteData.append([])
	arrWriteData[i].append(d)

with open('ber.csv', 'w') as f:
	wr = csv.writer(f)
	wr.writerows(arrWriteData)

BER(x)の出力です。1MUI(=1e6UI)のデータなのでy軸は1e-6で切れています。

f:id:kent_s:20180817172225p:plain

左側のバスタブを拡大。
f:id:kent_s:20180817172340p:plain

BER=1e-5の後半あたりから曲線が荒い。
これもサンプル数が1MUIの為、BER=5e-5付近で精度が落ちている為。
x=7psec(赤で示したあたり)がDjが無くなるポイントだが非常に滑らか。
つまり、このプロットからRjを予測するのは困難。

2.2 Qプロットの作成

2.1で作ったber.csvからQプロットを作成。
信号遷移率ρTは1とした。
今回のデータは毎サイクルにジッタが発生しているので、信号遷移率は100%となる。

import numpy as np
import csv 
import sys
import matplotlib.pyplot as plt
from scipy import special

fltUi = 100e-12
intRange = 1000
fltRho = 1

args = sys.argv
path = args[1]

#Read BER data
arrBer = []
with open(path, 'r') as f:
	rr = csv.reader(f)
	for rd in rr:
		arrBer.append(float(rd[0]))

#Calc Q plot
arrQ = []
for ber in arrBer:
	q = np.sqrt(2) * special.erfinv(1 - 1/fltRho*ber)
	arrQ.append(q)

#Generate Time axis
arrT = np.arange(0, fltUi, fltUi/intRange)

#Plot Q
plt.figure()
#plt.scatter(np.arange(intRange), arrQ)
plt.plot(arrT, arrQ)
plt.ylim(8, 0)
plt.show()

#Output Data
arrWriteData = []
for i,d in enumerate(arrQ):
	arrWriteData.append([])
	arrWriteData[i].append(d)

with open('q.csv', 'w') as f:
	wr = csv.writer(f)
	wr.writerows(arrWriteData)

Q(x)の出力。
f:id:kent_s:20180817172742p:plain

左のバスタブの拡大。
f:id:kent_s:20180817172849p:plain

バスタブカーブが線形に近いことが確認できる。この線形性を利用してσを計算できる。
x=7psecのSjが無くなるポイントだが、変曲点よりも外側に来ている。
これは参考資料にある通り、Djのpp点Dj(pp)よりも、正規分布中心Dj(δδ)の方が内側に来ることが理由。

2.3 σの計算

最後にσの計算を行う。
BERのバスタブ底は精度が無いので、BER=2e-5のところで近似することにした。
BER=2e-5 ⇒ Q=4.265 (信号遷移率=1の時)なので、Q=4.265のところでfittingを取る。

import numpy as np
import csv 
import sys
import matplotlib.pyplot as plt
from scipy import special

fltUi = 100e-12
intRange = 1000
fltRho = 1

args = sys.argv
path = args[1]

#Read Q data
arrQ = []
with open(path, 'r') as f:
	rr = csv.reader(f)
	for rd in rr:
		arrQ.append(float(rd[0]))

#Search Q 
fltQs = 4.265

def getNearestIndex(list, val):
	return (abs(np.asarray(list) - val)).argmin()

intIlt = getNearestIndex(arrQ[:intRange/2], fltQs)
intIrt = getNearestIndex(arrQ[intRange/2:], fltQs) + intRange/2

#Generate Time axis
arrT = np.arange(0, fltUi, fltUi/intRange)

#Fitting
FitLt = np.polyfit(arrT[intIlt-2:intIlt+2], arrQ[intIlt-2:intIlt+2], 1)
FitRt = np.polyfit(arrT[intIrt-2:intIrt+2], arrQ[intIrt-2:intIrt+2], 1)

print FitLt
print FitRt

#Plot Q
plt.figure()
plt.scatter(arrT, arrQ, c='red')
plt.plot(arrT, np.poly1d(FitLt)(arrT))
plt.plot(arrT, np.poly1d(FitRt)(arrT))
plt.ylim(8, 0)
plt.xlim(0, fltUi)
plt.show()

QプロットとQ=4.265での線形近似直線。
f:id:kent_s:20180818054121p:plain

線形近似の係数はax+bとすると、バスタブ左側はa=9.19e11, b=-5.35. 右側はa=-1.04e12, b=97.0となった。

1/aがRjのσとなるが、計算すると、σ@左=1.1ps, σ@右=1.0psと計算できた。
これは1.1でRjを生成した時の1ps RMSと一致しており、Dual Dirac Methodで正しく見積もれていることが確認できた。

Dj(δδ)=-b/aなので、Dj(δδ)@左=5.8ps, Dj(δδ)@右=-6.3psとなる。
ここからBER=1e-12として、Tj=Dj(δδ)+7*σで見積もると、それぞれ13ps, -13psで、トータル26psとなる。
Dj(pp)=7psだったので、Tj=Dj(pp)+7σで見積もった場合は±14psとなるが、それよりも1psずつ悲観性を無くすことができている。