フーリエ変換と線形システムの基礎 (2) - 例題

使ってみるのが理解の早道ですから、解説編の話題を手軽に体験できるような例題を用意してみました。

まず実験のための道具を用意します。

1. 準備

1.1. 離散フーリエ変換プログラム

まず、手軽に離散フーリエ変換が使えるように、fft という名前のプログラムを作っておきましょう。次のような仕様で設計します。

  コマンド	fft - FFT による離散フーリエ変換 (DFT)
  構文		fft [オプション] [入力ファイル [出力ファイル]]
  オプション	-B	入力データに Blackman-Harris フィルタを掛けてから変換
		-F	入力データに flat-top フィルタを掛けてから変換
		-M	入力データに Hamming フィルタを掛けてから変換
		-N	入力データに Hanning フィルタを掛けてから変換
		-R	入力データに Riesz フィルタを掛けてから変換
		-a	出力データを絶対値と偏角で表示
		-i	フーリエ逆変換(IDFT)を実行
		-p	出力データをパワー(dB)で表示
		-r	FFT 演算のみ実行
		-s	ステップ状関数補正を実行
		-x	-i を指定したとき、実部だけを表示

入力データは下記の形式とします。 実部虚部 # -r を指定したとき時刻実部 [虚部] # -i を指定しないとき (虚部 0 は省略可能) 周波数実部虚部 # -i を指定したとき -r を指定したときは FFT 演算のデータ順で全てのデータを入力しますが、-i を指定した場合は後半の負の周波数部分を省略することにして、fft プログラム内部で後半を埋めます。

出力データも同じ形式で、-r か -i を指定しない場合は、後半の負の周波数領域を出力しません。

自作が面倒なかたはfftプログラムをご利用ください。文字コードは EUC です。Chirp z-Transform 等、上記以外の機能も含んでいます。

1.2. 可視化ツール

波形やスペクトルを「見える化」するためのグラフ化ツールが必要ですが、良く使われる gnuplot は標準入力からデータが読めない点で Unix 的でなく、少し不便を感じますので、私自身は出力データをパイプでグラフ化ツールに繋げば適当なグラフを描いてくれるようなコマンドを自作して使ってきました。

ここでは、自作のグラフ化ツールと同じような操作で使える簡単なスクリプト「p」を用意して、それを使うことにします。次のような仕様です。

  コマンド	p - グラフ生成
  構文		p [オプション]
  オプション	-x	x 軸データなし

標準入力から読み取ったデータの最初のフィールドを x 軸の値、2 番目以降のフィールドを y 軸の値として折れ線グラフを生成します。-x を指定した場合は 0, 1, 2, .. の x 軸値が省略されたものとみなします。

2. FFT 演算と DFT/IDFT の関係

通常 FFT のプログラムでは

                  N-1   
  A(k) + j*B(k) = Σx(j)*exp(j*2*π*k/N)
                  j=0
  ここに、
	j = sqrt(-1)

を実行しますから、離散フーリエ変換、離散フーリエ逆変換の間には、

  F(k) = T*A(k) -j*T*B(k)
  x(k) = F*A(k) -j*F*B(k)
  ここに、
	T = 時間軸のサンプリング間隔
	F = 周波数軸のサンプリング間隔
	F(k) = x(k) の離散フーリエ変換
	x(k) = F(k) の離散フーリエ逆変換

の関係があることと、FFT 演算で得られた結果がメモリ上で下記のように並んでいることを理解しておく必要があります。

  実部 : R(0) R(1) R(2) ... R(N-1) R(N) R(-N+1) R(-N+2) .. R(-2) R(-1)
  虚部 : X(0) X(1) X(2) ... X(N-1) X(N) X(-N+1) X(-N+2) .. X(-2) X(-1)

中央の次からの後半は負の周波数のデータです。 x(j) が実時間関数なら、R(k) は偶関数、X(k) は奇関数ですから、

  R(k) = R(N - k)   (0 < k <= N/2)
  X(k) = -X(N - k)   (0 < k <= N/2)
  X(N/2) = 0

になります。

X(k) の刻みを半分にしたいときは、x(j) の後に 0 を N 個追加して、2*N 個のデータを FFT すればよいのですが、これが Bluestein's FFT algolithm (chirp z-transform) で使われています。

極めて有用で離散フーリエ変換には無難な減衰指数関数

  f(t) = exp(-t)
  F(f) = 1 / (1 + j*2*π*f)
  ここに、
	t = 時間 (s)
	f = 周波数 (Hz)

を使って、上記 FFT 演算の性格を確認してみましょう。

N=32	# サンプル数
T=0.25	# サンプリング間隔
awk '
BEGIN {
  N = '"$N"'
  T = '"$T"'
  for (i = 0; i < N; i++)
        print (i == 0) ? 0.5 : exp(-T * i), 0
  exit 0
}' | fft -r | awk '
BEGIN {
  N = '"$N"'
  T = '"$T"'
  F = 1 / N / T	# 周波数間隔
  PI2 = atan2(0, -1) * 2
}
{
  # FFT -> DFT
  $1 *= T
  $2 *= -T
  # 真値
  w = PI2 * (NR - 1) * F
  x = 1 / (1 + w * w)
  y = -w * x
  print $1, $2, $3, x, y
}
'

i = 0 の不連続点では (f(0) + f(∞)) / 2 にしています。真値である f(0) = 1 にした場合と比較してください。

逆変換の場合は、次のようになります。

N=32
F=0.125
awk '
BEGIN {
  N = '$N'
  F = '$F'
  PI2 = atan2(0, -1) * 2
  for (i = 0; i <= N / 2; i++) {
	w = PI2 * F * i
	x[i] = 1 / (1 + w * w)
	y[i] = -w * x[i]
  }
  y[N / 2] = 0	# N/2 で点対象にする
  for (i = 0; i <= N / 2; i++)
	print x[i], y[i]
  for (i = 1; i < N / 2; i++)
	print x[N / 2 - i], -y[N / 2 - i]
  exit 0
}' | fft -r | awk '
BEGIN {
  F = '$F'
  N = '$N'
  T = 1 / F / N
}
{
  # FFT -> IDFT
  $1 *= F
  $2 *= -F
  # 真値
  t = T * (NR -1)
  x = exp(-t)
  y = 0
  print $1, $2, x, y
}
'

スペクトルの後半は N/2 で折り返しますが、虚部については符号を反転し、点対象にするため N/2 の値を 0 にしています。

逆変換の終りのほうで振動を生じて真値との誤差が増えていますが、これが折返し歪 (alias) です。周波数間隔を減らしてデータ数を増やせば、さらに精度が上がります。

y[N/2] を 0 にしないと逆変換で虚部の誤差が増えることを確認してください。

FFT 演算と離散フーリエ変換/離散フーリエ逆変換の関係が理解できたら、以下、スケーリングや虚部の符号反転、時間軸や周波数軸データの生成は fft プログラムに任せることにして先に進みます。

3. 離散フーリエ変換

3.1. デルタ関数

N=32
/usr/bin/awk '
BEGIN {
  N = '"$N"'
  for (i = 0; i < 1; i += 1 / N)
        if (i == 0)
                print i, N      # 面積
        else
                print i, 0
  exit 0
}' | fft | awk '
BEGIN {
  N = '"$N"'
}
{
  x = 1
  y = 0
  print $1, $2, $3, x, y
}
'

f(0) は変換区間の面積にしておかないとスペクトルのスケールが変わってしまいます。

この結果を「| fft -i」に接続すると、元に戻ります。

3.2. 定数

N=32
/usr/bin/awk '
BEGIN {
  N = '"$N"'
  for (i = 0; i < N; i++) {
        t = i / N
        print t, 1
  }
  exit 0
}' | fft

この結果を「| fft -i」に接続すると、元に戻ります。

周波数 0 の成分しかないのは当然ですが、波高値を変えてみてください。

また、ステップ関数補正を行うと、ステップ関数のスペクトルの近似値が得られます。

3.3. 正弦波

# cos(t)
/usr/bin/awk '
BEGIN {
  pi = atan2(0, -1)
  N = 32
  for (i = 0; i < 1; i += 1 / N)
	print i, cos(2 * pi * i)
  exit 0
}' | fft

# sin(t)
/usr/bin/awk '
BEGIN {
  pi = atan2(0, -1)
  N = 32
  for (i = 0; i < 1; i += 1 / N)
	print i, sin(2 * pi * i), 0
  exit 0
}' | fft

いずれも、「| fft -i」で元に戻りますが、例えば、

  print i, sin(2.3 * pi * i), 0

等、周期に合わせて波形を切り取らないと、折返し歪が発生することを確認しましょう。

3.4. 指数減衰関数

C-R 直列回路にステップ電圧を加えたときの R の電圧をフーリエ変換して、真の値と比較しています。このような減衰波形だと小さな誤差で済みます。不連続点 t = 0 の値を 1 でなく 0.5 にすれば、t が小さいときの誤差が大幅に減ります。

R=100
C=50e-12
N=32
T=1e-9
awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
  for (i = 0; i < N; i++) {
	t = i * T
	f = exp(-t / R / C)
	if (i == 0) f = 0.5	# 不連続点 - t が小さいときの誤差が減る
	print t, f, 0
  }
  exit 0
}' | fft | awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
  PI2 = atan2(0, -1) * 2
}
{
  # 真値
  w = PI2 * $1 * C * R
  x = R * C / (1 + w * w)
  y = -w * x
  print $1, $2, $3, x, y
}
'

「| fft -i」に繋いで逆変換してみましょう。

4. ステップ状波形

0 から始まり 0 に戻らない関数を離散フーリエ変換すると、0 に戻る周期関数として扱われ、大きな誤差を生じますが、ある時刻以降の値が一定なら、ステップ状関数 (Step-like function) 補正が使えます。

4.1. ステップ関数

N=32
/usr/bin/awk '
BEGIN {
  N = '"$N"'
  for (i = 0; i < N; i++) {
	t = i / N
	print t, 1
  }
  exit 0
}' | fft -s | awk '
BEGIN {
  N = '"$N"'
  PI2 = atan2(0, -1) * 2
  t0 = 1 / N
}
{
  # 真値
  w = PI2 * $1
  if (w == 0) {
	x = 1	# 面積
	y = 0	# 本当は負の無限大だが、ステップ補正は 0 を扱わない
  }
  else {
	x =  0.5 / N
	y = -1 / w
  }
  print $1, $2, $3, x, y
}
'

ステップ補正を行わないと、虚部のスペクトルがまるで変わってしまいます。

4.2. 指数増加関数

C-R 直列回路にステップ電圧を加えたときの C の電圧をフーリエ変換して、真の値と比較しています。ステップ状関数補正の有効性がよくわかります。

R=100
C=50e-12
N=32
T=1e-9
awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
  for (i = 0; i < N; i++) {
	t = i * T
	f = 1 - exp(-t / R / C)
	print t, f, 0
  }
  exit 0
}' | fft -s | awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
  PI2 = atan2(0, -1) * 2
}
{
  w = PI2 * $1
  if (w == 0)
	next
  x = w * C * R
  y = -1 / w / (1 + x * x)
  x *= y
  print $1, $2, $3, x, y
}
'

4.3. 傾斜関数

傾斜関数のフーリエ変換結果を真値と比較しています。「fft -s 」を「fft」にしてステップ状関数補正をやめると、大きな誤差が出てきます。

N=256
T=1
T0=0.2
T0=0.1
awk '
BEGIN {
  T = "'"$T"'" + 0
  t0 = "'"$T0"'" + 0
  N = "'"$N"'" + 0
  dt = T / N
  for (i = 0; i < N; i++) {
	t = i * dt
	if (t < t0)
		f = t / t0
	else
		f = 1
	print t, f
  }
  exit 0
}
' | fft -s | awk '
BEGIN {
  N = '"$N"'
  t0 = '"$T0"'
  PI2 = atan2(0, -1) * 2
}
{
  # 真値
  w = PI2 * $1
  if (w == 0) {
	x = N * (1 - t0 / 2) / N
	y = 0
  }
  else {
	a = 2 * sin(w * t0 / 2) / t0 / w / w
	x = -a * sin(t0 / 2 * w)
	y = -a * cos(t0 / 2 * w)
  }
  print $1, $2, $3, x, y
}
'

4.4. 少し複雑なステップ状波形

下記の論文で使われた面白い波形です。実数だけで計算するのは面倒ですが、複素数で扱えば簡単です。

  G.D.Cormac and J.N.McMullin,- Frequency De-Aliased FFT Analysis of Step-Like Functions
  IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT. VOL.40, NO. 4. AUGUST 1991

どんな波形かは fft コマンドに入れる前のデータをグラフ化してみてください。

# f(t) = { Us + U0*exp(-a(t - ts))*sin(2*π*fx*(t - ts)))*u(t - ts)
#	+ 0.5*Us*(1 - cos(π*t/ts))*(u(t) - u(t - ts))
# ts = 200*T/N, N = 512, a = 12.5/(N - 1), fx = 75 / (N - 1)
# Us = 1, U0 = 0.125
N=512
Us=1
U0=0.125
A=12.5
T=1
S=0.4
F=75
/usr/bin/awk '
BEGIN {
  N = "'"$N"'" + 0
  Us = "'"$Us"'" + 0
  U0 = "'"$U0"'" + 0
  A = "'"$A"'" + 0
  F = "'"$F"'" + 0
  T = "'"$T"'" + 0
  S = "'"$S"'" + 0
  a = A
  fx = F
  ts = S * T
  dt = T / N
  PI = atan2(0, -1)
  for (i = 0; i < N; i++) {
	t = i * dt
	if (t <= ts) {
		f = 0.5*Us*(1 - cos(PI*t/ts))
	}
	else {
		f = Us + U0*exp(-a*(t - ts))*sin(2*PI*fx*(t - ts))
	}
	print t, f
  } 
  exit 0
}
' | fft -s | /usr/bin/awk '
BEGIN {
  N = "'"$N"'" + 0
  Us = "'"$Us"'" + 0
  U0 = "'"$U0"'" + 0
  A = "'"$A"'" + 0
  F = "'"$F"'" + 0
  T = "'"$T"'" + 0
  S = "'"$S"'" + 0
  a = A
  fx = F
  ts = S * T
  dt = T / N
  df = 1 / N / dt
  PI = atan2(0, -1)
  PI2 = atan2(0, -1) * 2
  k = 10 / log(10)
}
{
  if (NR == 1)
	next
  true = True()
  dealias = Dealias()
  print $1, dB($2, $3), true, dealias
}
# F(ω) = Us*exp(-j*ω*ts)/(j*ω)
#	+ U0*2*π*fs*exp(-j*ω*ts)/((a + j*ω)^2 + (2*π*fx)^2))
#	+ 0.5*Us*(1 - exp(-j*ω*ts)/(j*ω)
#	+ 0.5*Us*j*ω*(1 + exp(-j*ω*ts))/(ω^2 - (π/ts)^2)
function True() {
  w = PI2 * $1
  Exp(0, -w*ts)
  u1 = x; v1 = y
  Div(Us*u1, Us*v1, 0, w)
  x1 = x; y1 = y;
  w1 = U0*PI2*fx
  Div(w1*u1, w1*v1, a*a - w*w + PI2*fx*PI2*fx, 2*a*w)
  x2 = x; y2 = y
  Div(0.5*Us*(1 - u1), -0.5*Us*v1, 0, w)
  x3 = x; y3 = y
  Mul(0, 0.5*Us*w, 1 + u1, v1)
  Div(x, y, w*w - PI/ts*PI/ts, 0)
  x4 = x; y4 = y
  x = x1 + x2 + x3 + x4; y = y1 + y2 + y3 + y4
  return dB(x, y)
}
# Fdealias(ω) = (sin(n/N)/(n/N))^2*F(ω)
function Dealias(  p, q) {
  p = (NR - 1) / N
  q = sin(p)
  p = (p == 0) ? 1 : q * q / (p * p)
  return dB(p * x, p * y)
}
function dB(x, y) {
  return k * log(x * x + y * y)
}
function Mul(x1, y1, x2, y2) {
  x = x1 * x2 - y1 * y2
  y = x1 * y2 + y1 * x2
}
function Div(x1, y1, x2, y2,  r2) {
  r2 = x2 * x2 + y2 * y2
  x = (x1 * x2 + y1 * y2) / r2
  y = (y1 * x2 - x1 * y2) / r2
}
function Exp(x1, y1,  a) {
  a = exp(x1)
  x = a * cos(y1)
  y = a * sin(y1)
}
'

5. 窓関数

2 つスペクトル成分を含んでいますが、近い場所なので、始端と終端の不連続から来る折返し歪に埋まって、パワースペクトルが分離できません。

しかし、Blackman-Harris 窓とか flat-top 窓を使うと分離できて、後者では成分の大きさの割合も正確にわかります。いくつかの窓関数を比較してみましょう。

T = 1 にすると、flat-top 窓でのみ分離できます。

N=256
awk '
BEGIN {
  N = '$N'
  T = 1.1
  a = 0.01
  PI2 = atan2(0, -1) * 2
  for (i = 0; i < N; i++) {
        t = i / N / T
        f = cos(PI2 * 5 * t) + a * sin(PI2 * 6 * t)
        print t, f
  }
  exit 0
}' | fft -pN

6. 離散フーリエ逆変換

6.1 指数減衰関数

C-R 直列回路にステップ電圧を印加したときの R の電圧を周波数領域で考えると、伝達関数(インパルス応答)を 1/(j*ω) 倍(積分))して、

  j*ω*R*C / (1 + j*ω*R*C) / (j*ω)

になりますから、これを離散フーリエ逆変換して、真値と比較します。

R=100
C=50e-12
N=256
T=1e-10
awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
  F = 1 / N / T
  PI2 = atan2(0, -1) * 2
  for (i = 0; i < N; i++) {
	f = F * i
	w = PI2 * f * C * R
	x = R * C / (1 + w * w)
	y = -w * x
	print f, x, y
  }
  exit 0
}
' | fft -i | awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
}
{
  # 真値
  t = $1
  f = exp(-t / R / C)
  print $1, $2, f
}
'

先端と最後で僅かな折り返し歪が出ますが、無難な関数ですから、うまくゆきます。

6.2 指数増加関数

R-C 直列回路にステップ電圧を印加したときの C の電圧を周波数領域で考えると、関数(インパルス応答)を 1/(j*ω) 倍(積分))して、

  1 / (1 + j*ω*R*C) / (j*ω)

になりますが、これを離散フーリエ逆変換すると出発点に戻る周期関数にされてしまいますから、変換区間外の部分に終点の波高値と同じ電圧の負のステップ関数を追加して、そのスペクトルを加算することで、出発点に戻らないステップ状波形の周期関数に変換してから離散フーリエ逆変換します。

R=100
C=50e-12
N=32
T=1e-9
awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
  t0 = T * (N - 1)	# 変換区間の終端
  F = 1 / N / T
  PI2 = atan2(0, -1) * 2
  for (i = 0; i < N; i++) {
	f = F * i
	w = PI2 * f
	if (w == 0) {
		x = t0 - C * R * (1 - exp(-t0 / C / R))	# 面積
		y = 0
	}
	else {
		x = w * C * R
		y = -1 / w / (1 + x * x)
		x *= y
		# t0 から立ち上がる負のステップ波形のスペクトルを加える
		x += sin(w * t0) / w
		y += cos(w * t0) / w
	}
	print f, x, y
  }
  exit 0
}
' | fft -i | awk '
BEGIN {
  R = '"$R"'
  C = '"$C"'
  N = '"$N"'
  T = '"$T"'
}
{
  # 真値
  t = $1
  f = 1 - exp(-t / R / C)
  print $1, $2, f
}
'

6.4 ステップ関数

ステップ関数を t0 だけ遅らせた周波数領域の関数

  F(f) = exp(-j*2*π*t0)/(j*2*π)

も前記の方法で扱うことができますが、ここでは

  F(k) = F(f - j*a/(2*π))/T  (k = 0, 1, 2, .. N/2)

の数列を離散フーリエ逆変換し、最後に exp(a*k) 倍する方法を使っています。

後半で大きく拡大されるため振動が出ますが、 Hanning 窓をかけてから離散フーリエ逆変換すると抑えられることを確認してください。

# Step 関数のフーリエ逆変換
# F(ω) = exp(-j*ω*t0)/(j*ω)
N=64
T=6.4
T0=1
A=6	# damping factor
/usr/bin/awk '
BEGIN {
  N = "'"$N"'" + 0
  T = "'"$T"'" + 0
  t0 = "'"$T0"'" + 0
  a = "'"$A"'" / T
  PI2 = atan2(0, -1) * 2
  for (t = 0; t <= N / 2; t++) {
	w = PI2 * t / T
	# exp(-t0 * (a + j * ω))
	e = exp(-a * t0)
	x1 = e * cos(w * t0)
	y1 = -e * sin(w * t0)
	# exp(-t0*(a+j*ω)/(a + j*ω) / T
	d = a * a + w * w
	x = (a * x1 + w * y1) / d / T
	y = (a * y1 - w * x1) / d / T
        print t, x, y
  }
  exit 0
}
' | fft -iN | /usr/bin/awk '
BEGIN {
  N = "'"$N"'" + 0
  T = "'"$T"'" + 0
  T0 = "'"$T0"'" * T
  a = "'"$A"'"
  PI2 = atan2(0, -1) * 2
  df = 1 / N / T
}
{
  t = $1 * T
  y = exp(a * $1) * $2
  print t, y
}
'

7. ラプラス逆変換

ステップ関数を t0 だけ遅らせたラプラス変換の像関数

  F(f) = exp(-t0*s)/s

を逆変換して時間関数に戻します。

N=64
T=6.4
T0=1
A=5	# damping factor
/usr/bin/awk '
BEGIN {
  N = "'"$N"'" + 0
  T = "'"$T"'" + 0
  t0 = "'"$T0"'" + 0
  a = "'"$A"'" / T
  PI2 = atan2(0, -1) * 2
  for (k = 0; k <= N / 2; k++) {
	f = k / T
	x1 = a; y1 = PI2 * f
	# exp(-(a + j*2*π*k/T)*t0)
	Exp(-x1 * t0, -y1 * t0)
	Div(x, y, x1, y1)
        print f, x, y
  }
  exit 0
}
function Div(x1, y1, x2, y2,  r2) {
  r2 = x2 * x2 + y2 * y2
  x = (x1 * x2 + y1 * y2) / r2
  y = (y1 * x2 - x1 * y2) / r2
}
function Exp(x1, y1,  a) {
  a = exp(x1)
  x = a * cos(y1)
  y = a * sin(y1)
}
' | fft -iN | /usr/bin/awk '
BEGIN {
  N = "'"$N"'" + 0
  T = "'"$T"'" + 0
  T0 = "'"$T0"'" * T
  a = "'"$A"'"
}
{
  # exp(a*k) 倍する
  y = $2 * exp(a / N * (NR - 1))
  print (NR - 1) * T / N, y
}
'

Hanning フィルタを外すとかなりの振動が出ます。

8. VNA (vector Network Analyzer) シミュレーション

8.1. S11 - 無損失ケーブル

75 Ohm の無損失同軸ケーブルを 50 Ohm の VNA (Vector Network Analyzer) に接続して S11 を測定するシミュレーションです。ケーブルの長さは 16 回分の反射波が見られるように選んでいます。

T=20e-12	# サンプリング間隔
N=512
/usr/bin/awk '
BEGIN {
  T = '$T'
  N = '$N'
  M = N / 16	# 16 回の反射まで観測する
  Z0 = 75	# ケーブルの特性インピーダンス
  r = (50 - Z0) / (50 + Z0)	# ケーブル端の反射係数
  a = b = Z0 / (50 + Z0)	# a = 始端電圧 b = 進行波
  c = 0.5	# 無反射始端電圧
  for (i = 0; i < N; i++) {
	if (i == 0)
		v = 0.5 * (a - c) / c		# 不連続点
	else if (i % M == 0) {
		a += b * r * (1 + r)
		b = b * r * r
		v =  (a - c) / c
	}
	else
		v =  (a - c) / c
	print T * i, v, 0
	v0 = v
  }
  exit 0
}' | fft | awk '
BEGIN {	# S11
  T = '$T'
  N = '$N'
  PI2 = atan2(0, -1) * 2
  k = 10 / log(10)
}
{
  if ($1 == 0)
	next
  # j*ω 倍(微分)してインパルス応答に変換
  x = -PI2 * $1 * $3
  y = PI2 * $1 * $2
  p = k * log(x * x + y * y)	# dB
#  print $1, p, atan2(y, x)
#  print $1, x, y
  print $1, p
}
'

8.2. S21 - インダクタ

VNA の Port1 と Port2 の間にインダクタを挿入した場合です。TDT 波形は指数関数的に増加しますから、観測を打ち切った後も増加し続けるため、ステップ状関数補正を行っても、後半の周波数ではかなりの誤差を生じます。

T=20e-12	# サンプリング間隔
N=512		# サンプリング数
L=4e-8		# インダクタンス (H)
R=50		# VNA のシステムインピーダンス
/usr/bin/awk '
BEGIN {
  T = '$T'
  N = '$N'
  L = '$L'
  R = '$R'
  for (i = 0; i < N; i++) {
	t = T * i
	print t, 1 - exp(-2 * R / L * t)
  }
  exit 0
}' | fft -s | awk '
BEGIN {	# S21
  T = '$T'
  N = '$N'
  L = '$L'
  R = '$R'
  PI2 = atan2(0, -1) * 2
  k = 10 / log(10)
}
{
  if ($1 == 0)
	next
  # j*ω 倍(微分)してインパルス応答に変換
  w = PI2 * $1
  x = -w * $3
  y = w * $2
  p = k * log(x * x + y * y)	# dB
  # 真値
  u = 2 * R
  v = w * L
  d = 2 / (u * u + v * v)
  X = R * u * d
  Y = -R * v * d
  q = k * log(X * X + Y * Y)	# dB
  print $1, p, q
#  print $1, p, atan2(y, x)
#  print $1, x, y
}
'

9. 注

9.1. グラフ化ツール

標準入力から読み込んだデータを gnuplot でグラフ化する簡単なスクリプト例です。

$ cat p
#! /bin/sh

while [ $# -gt 0 ]
do
  if [ $1 = "-L" -o $1 = "-Lxy" ]
  then
	SCALE="set logscale"
  elif [ $1 = "-Lx" ]
  then
	SCALE="set logscale x"
  elif [ $1 = "-Ly" ]
  then
	SCALE="set logscale y"
  elif [ $1 = "-Tpng" ]
  then
	TERMINAL="set terminal pngcairo"
  fi
  shift
done

/usr/bin/awk '
BEGIN {
  SCALE = "'"$SCALE"'"
  TERMINAL = "'"$TERMINAL"'"
  nf = N = L = 0
  if (TERMINAL != "") {
	print TERMINAL
  }
  print "set multiplot"
  if (SCALE != "")
	print SCALE
}
/^#/ {
  next
}
/^$/ {
  if (N > 0)
	plot()
  next
}
/^s/ {
  print "set title \047" substr($0, 3) "\047"
  next
}
/^v/ {
  print "set origin", $2 "," $4
  print "set size", ($3 - $2) "," ($5 - $4)
  next
}
/^wx/ {
  if (NF == 3) {
	print "set logscale x"
	print "set xtics 10"
	print "set xrange [" 10^$2 ":" 10^$3 "]"
	next
  }
  print "set xrange [" $2 ":" $3 "]"
  if (NF == 5) {
	print "set xtics", $4
	print "set mxtics", $4 / $5
  }
  else if (NF == 4)
	print "set xtics", $4
  next
}
/^wy/ {
  if (NF == 3) {
	print "set logscale y"
	print "set yrange [" 10^$2 ":" 10^$3 "]"
	print "set ytics 10"
	next
  }
  print "set yrange [" $2 ":" $3 "]"
  if (NF == 5) {
	print "set ytics", $4
	print "set mytics", $4 / $5
  }
  else if (NF == 4)
	print "set ytics", $4
  next
}
/^x/ {
  print "set xlabel \047" substr($0, 3) "\047"
  next
}
/^y/ {
  print "set ylabel \047" substr($0, 3) "\047"
  next
}
/^l/ {
  label[L++] = substr($0, 3)
  next

}
{
  for (i = 1; i <= NF; i++)
	d[N, i] = $i
  N++
  if (nf == 0)
	nf = NF
  else if (nf != NF)
	print "abnormal record", $0 >"/dev/stderr"
}
END {
  if (N > 0)
	plot()
}
function plot(  i, j, k) {
  k = 0
  printf "plot \"-\" title \"" label[k++] "\" with lines"
  w = ""
  for (i = 2; i < nf; i++)
	printf ", \"-\" title \"" label[k++] "\" with lines"
  printf "\n"
  for (i = 2; i <= nf; i++) {
	for (j = 0; j < N; j++)
		print d[j, 1], d[j, i]
	print "e"
  }
  for (i = 0; i < L; i++)
	label[i] = ""
  w = ""
  nf = N = L = 0
  print "reset"
}
' | gnuplot -persist

使いかたは、例えば、

  {fft 演算を使う実験用スクリプト} | p

でグラフが出てきますので、見終ったらグラフ画面にカーソルを移動して「q」キーを押して終了します。

9.2. 複素関数

フーリエ変換は複素関数ですから、式の計算も複素数でやるほうが楽です。c++ ならいとも簡単ですが、awk を使うときは戻り値に x, y の広域変数を使うように設計するのが良いと思います。

手間を省くため、よく必要になる複素関数をまとめておきます。

# complex functions
function Add(x1, y1, x2, y2) {
  x = x1 + x2
  y = y1 + y2
}
function Sub(x1, y1, x2, y2) {
  x = x1 - x2
  y = y1 - y2
}
function Mul(x1, y1, x2, y2) {
  x = x1 * x2 - y1 * y2
  y = x1 * y2 + y1 * x2
}
function Div(x1, y1, x2, y2,  r2) {
  r2 = x2 * x2 + y2 * y2
  x = (x1 * x2 + y1 * y2) / r2
  y = (y1 * x2 - x1 * y2) / r2
}
function Exp(x1, y1,  a) {
  a = exp(x1)
  x = a * cos(y1)
  y = a * sin(y1)
}
function Log(x1, y1) {
  x = 0.5 * log(x1 * x1 + y1 * y1)
  y = atan2(y1, x1)
}
function Sqrt(x1, y1,  a) {
  x = 0.5 * log(x1 * x1 + y1 * y1) * 0.5
  y = atan2(y1, x1) * 0.5
  a = exp(x)
  x = a * cos(y)
  y = a * sin(y)
}
function Pow(x1, y1, x2, y2,  x3, y3) {
  Log(x1, y1)
  x3 = x
  y3 = y
  Mul(x3, y3, x2, y2)
  Exp(x, y)
}
function Sin(x1, y1,  e, f) {
  e = exp(y1)
  f = 1 / e
  x = 0.5 * cos(x1) * (e - f)
  y = 0.5 * sin(x1) * (e + f)
}
function Cos(x1, y1,  e, f) {
  e = exp(y1)
  f = 1 / e
  x = 0.5 * cos(x1) * (f + e)
  y = 0.5 * sin(x1) * (f - e)
}
function Tan(x1, y1,  e, f, d) {
  e = exp(2 * y1)
  f = 1 / e
  d = cos(2 * x1) + 0.5 * (e + f)
  x = sin(2 * x1) / d
  y = 0.5 * (e - f) / d
}
function Sinh(x1, y1,  e, f) {
  e = exp(x1)
  f = 1 / e
  x = 0.5 * (e - f) * cos(y1)
  y = 0.5 * (e + f) * sin(y1)
}
function Cosh(x1, y1,  e, f) {
  e = exp(x1)
  f = 1 / e
  x = 0.5 * (e + f) * cos(y1)
  y = 0.5 * (e - f) * sin(y1)
}
function Tanh(x1, y1,  e, f, d) {
  e = exp(2 * x1)
  f = 1 / e
  d = 0.5 * (e + f) + cos(2 * y1)
  if (d == 0) {
	e = exp(x1)
	f = 1 / e
	x = 0
	y = (e + f) * sin(y1)
  }
  else {
	x = 0.5 * (e - f) / d
	y = sin(2 * y1) / d
  }
}