フーリエ変換と線形システムの基礎 (1) - 解説

フーリエ変換(Fourier Transform)は理工学の多くの分野で重要な役割を果たしますが、 ここでは主として工学分野で線形システムの解析に必要となる、時間領域と周波数領域 の変換を中心に、実務的かつ基 礎的事項を解説します。実例については続編を見て ください。

フーリエ変換で行き来できる時間領域と周波数領域は同じ実体の両面ですから、同じ情報 を含んでいるにもかかわらず、目に見えるものはまったく 別で、それぞれの長所があります。 つまり、トヨタ流発想から見れば、 フーリエ変換は「見える化」 の重要な手段の1つなのです。

1. (連続)フーリエ変換

まず、(連続)フーリエ変換の要点を解説しますが、 詳細は、下記の文献などを参照してください。

  Athanasios Papoulis,- THE FOURIER INTEGRAL AND ITS APPLICATIONS
	(McGraw-Hill)
何ともすばらしい本で、 「工学のための 応用フーリエ積分」(オーム社)という和訳もあります。 また、 http://www.ingelec.uns.edu.ar/pds2803/Materiales/LibrosPDF/Papoulis/TOC.HTM でも見られます。

1.1. フーリエ変換の定義

時間領域の関数 f(t) と周波数領域の関数 F(f) を相互に変換する フーリエ変換の定義にはいくつかの流儀がありますが、 ここでは、角周波数(ω)でなく周波数(f)を使う、下記の形式を使うことにします。 (注1)

          ∞
  F(f) = ∫f(t)*exp(-j*2*π*f*t)*dt                        (1.1)
        -∞
       = R(f) + j*X(f)
       = A(f)*exp(j*φ(f))
          ∞
  f(t) = ∫F(f)*exp(j*2*π*f*t)*df                         (1.2)
        -∞
  ここに、
	t = 時間 (s)
	f = 周波数 (Hz)
	f(t) = 時間を変数とする関数
	F(f) = f(t) のフーリエ変換
	R(f) = F(f) の実部
	X(f) = F(f) の虚部
	A(f) = sqrt(R(f)^2+X(f)^2)
	φ(f) = atan(X(f)/R(f))
	j = sqrt(-1)

f(t) を無数の指数関数の和に分解した結果が F(f) で、 F(f) を f(t) のフーリエスペクトル (Fourier Spectrum) と呼ぶのが普通です。

1.2. フーリエ変換対の基礎的な性質

以下、f(t) と F(f) の対応を

  f(t) <-> F(f)
と書くことにすれば、下記の関係は応用範囲が広く、重要です。なかでも、 畳込みは図でも書きながらゆっくり考えないと理解しにくいのですが、 特別重要なものの1つで、理解の労力の投資のしがいがあります。
  1. 線形性 - 線形システム解析の基盤
      a1*f1(t) + a2*f2(t) <-> a1*F1(f) + a2*F2(f)   (1.2.1)
    
  2. 対称性
      f(t) <-> F(-f)    (1.2.2)
    
  3. 時間軸の伸縮 (abs(x) は x の絶対値)
      f(a*t) <-> (1/abs(a))*F(f/a)  (1.2.3)
    
  4. 周波数軸の伸縮
      (1/abs(a))*f(t/a)h(a*t) <-> F(a*f)    (1.2.4)
    
  5. 時間軸の推移
      f(t - t0) <-> F(f)*exp(-j*2*π*f*t0)    (1.2.5)
    
  6. 周波数軸の推移 - 変調 (ラプラス逆変換でも重要)
      f(t)*exp(j*2*π*t*f0) <-> F(f - f0)    (1.2.6)
    
  7. 畳込み(重畳積分)定理 (convolution theorem) - h(t) がインパルス関数の場合は特に重要
            ∞
      y(t) = ∫x(t)*h(t - τ)*dτ <-> X(f)*H(f)    (1.2.7)
          -∞
      ここに、x(t) <-> X(f),  h(t) <-> H(f)
    
  8. 時間微分 - 微分演算の代数化
      (d/dt)^n(f(t)) <-> (j*2*π*f)^n*F(f)    (1.2.8)
    
  9. 周波数微分
      (-j*2*π*t)^n*f(t) <-> (d/df)^n(F(f))    (1.2.9)
    
  10. 時間積分 - 積分演算の代数化
       ∫∫..∫f(t)(dt)^n <-> F(f)/(j*2*π*f)^n    (1.2.10)
    
  11. 周波数積分
       (-j*2*π*t)^-n*f(t)(dt)^n <-> ∫∫..∫F(f)(df)^n    (1.2.11)
    
  12. 共役関数 (<f(x)> は f(x) の共役関数)
      <f(t)> <-> <F(-f)>     (1.2.12)
    
  13. モーメント定理 - M0 は面積
      (-2*π*j)^n*Mn = (d/df)^n(F(0))  (n = 0, 1, 2, ..)    (1.2.13)
      ここに、
            ∞
      Mn = ∫t^n*f(t)*dt  (n = 0, 1, 2, ..)
          -∞
    
            ∞
      M0 = ∫f(t)*dt = F(0)    (1.2.14)
          -∞
    
  14. エネルギ定理 - エネルギ保存
       ∞            ∞
      ∫f(t)^2*dt = ∫F(f)^2*df    (1.2.15)
      -∞          -∞
    

1.3. フーリエ変換対の例

以下の解説で必要になるフーリエ変換対をまとめておきます。 これだけでも、前記の基本的性質と組み合わせて、 かなりのことができます。 ステップ関数を微分(j*2*π*f倍)すると、 インパルス関数になることに注意してください。

  1. インパルス関数
      δ(t) <-> 1
      δ(t - t0) <-> exp(-j*2*π*f*t0)   (時間軸推移)
      1 <-> δ(f)
      ここに、
    	  ∞
    	 ∫δ(t)*dt = 1     δ(t) = 0 (t != 0)
    	-∞
      なお、
    	f(t)*δ(t) = f(0)*δ(t)
    	t*δ(t) = 0
    
  2. ステップ関数
      U(t) <-> δ(f)/2 + 1/(j*2*π*f)
      ここに、
    	U(t) = 1 (0 < t)
    	     = 0 (t < 0)
      なお、
    	(d/dt)U(t) = δ(t)
    
  3. 指数減衰関数
      exp(-a*t)*U(t) <-> 1/(a + j*2*π*f)
    
  4. 正弦波
      sin(2*π*f0*t) <-> (j/2)*(δ(2*π*f + f0) - δ(2*π*f - f0))
      cos(2*π*f0*t) <-> (j/2)*(δ(2*π*f + f0) + δ(2*π*f - f0))
    
  5. 傾斜関数
      V(t) <-> δ(f)/2 + sin(π*f*T)/(j*4*π*T*(π*f)^2)
      ここに、
    	V(t) = 0 (t < -T/2)
    	     = (t + T/2)/T (-T/2 <= t <= T/2)
    	     = 1 (T/2 < t)
    

2. 線形システム

線形時間的に不変 なシステムに f(t) を入力したとき、 出力 g(t) は下記のように表すことができます。(注2)

          ∞                     ∞
  g(t) = ∫f(t)*h(t - τ)*dτ = ∫f(t - τ)*h(τ)*dτ     (2.1)
        -∞                    -∞
  ここに、
	h(t) = インパルス関数
つまり、 h(t) だけでシステムのすべてがわかるという驚くべき 結論ですが、 これが線形回路理論の基盤になります。

g(t) <-> G(f), h(t) <-> H(f) とすれば、畳込み定理を使って、 周波数領域でシステムを決定することができます。

  G(f) = F(f)*H(f)                                        (2.2)
  ここに、
	H(f) = システム関数
	     = R(f) + j*X(f)
	     = A(f)*exp(j*φ(f))
	α(f) = -log(A(f)) = 減衰量 (neper)
	β(f) = -φ(f) = 位相遅れ (rad)
ここで、f(t) をインパルス関数にすれば、F(f) = 1 ですから、H(f) が がシステムのインパルス応答のフーリエ変換になることがわかります。

線形で時間的に不変なシステムの、もう1つの重要な性格は、 入力が指数関数 exp(j*ω*t) なら、 その出力もまた入力に比例した指数関数になる(注3)ことで、 f(t) と g(t) が指数関数の場合、H(f) の逆数は 「インピーダンス」(impedance) と呼ばれ、 電気工学や機械工学では極めて有効な概念になります。

例えば、R-C 直列回路にステップ電圧を印加したときの応答(C の電圧)(注4)は、

  f(t) = 1 - exp(t/(C*R))  (0 <= t)
  ここに、
	f(t) = R-C 直列回路に 1 V のステップ電圧を印加したときの C の電圧 (V)
	C = キャパシタンス (F)
	R = 抵抗 (Ohm)
	t = 時間 (s)
ですが、右辺をフーリエ変換すると、
  F(f) = δ(f)/2 + 1/(j*2*π*f*(1 + j*2*π*f*R*C))
ステップ関数のフーリエ変換を F(f) で割ればシステム関数が得られますが、結果は、
  H(f) = 1/(1 + j*2*π*f*R*C)
で、周波数領域の交流理論による計算と一致します。

なお、ここで、 ステップ応答を j*2*π*f 倍するとインパルス応答が得られたわけですが、 この問題は後で TDR/TDT 波形を扱うとき再考します。

3. 離散フーリエ変換と FFT

数値データを扱う場合はフーリエ変換の離散化が必要ですが、 FFT (Fast Fourie Transform) と呼ばれる優れたアルゴリズムが開発され、 効率のよい計算ができます。 この分野の詳細は、例えば、下記をご覧ください。

  E.Oran Brigham,- THE FIRST FOURIER TRANSFORM
	(Prentice-Hall)
明解な解説で、「高速フーリエ変換」(科学技術出版社)の和訳もあります。

3.1. 離散フーリエ変換による連続フーリエ変換の近似

(連続)フーリエ変換とフーリエ逆変換の対、

          ∞
  F(f) = ∫f(t)*exp(-j*2*π*f*t)*dt
        -∞
          ∞
  f(t) = ∫H(f)*exp(j*2*π*f*t)*dt
        -∞
を、下記の離散フーリエ変換と離散フーリエ逆変換の対で近似することを考えます。
                N-1
  F(n*dF) = dT*(Σ(f(k*dT)))*exp(-j*2*π*n*k/N)  (n = 0, 1, .., N - 1)  (3.1.1)
                k=0

                N-1
  f(k*dT) = dF*(Σ(F(n*dF)))*exp(j*2*π*n*k/N)  (k = 0, 1, .., N - 1)   (3.1.2)
                n=0

  ここに、
	N = 標本数
	dT = 時間軸上の標本化間隔
	dF = 周波数軸上の標本化間隔 = 1/(N*dT)
	T = N*dF = 切取り時間幅

離散化で不可避の 標本化誤差に加えて、 無限の時間を対象とする連続フーリエ変換と 有限な時間を対象とする離散フーリエ変換には大きな違いがあって、 下記の点が重要です。

  1. F(f) = 0 (fc <= f) なら、f(t) は dT = 1/(2*fc) の標本値だけから決定できる。 (Sampling Theorem - 標本化定理)
    サンプリング間隔が 1/(2*fc) を越えると、 F(f) に折返し歪(aliasing distortion)を生じ、 正しい結果が得られなくなります。
  2. 離散フーリエ変換は時間領域と周波数領域に於ける 周期性が前提で、 時間領域に於ける N 個の標本が周期関数の 1 周期になる。
    非周期関数の F(f) は f(t) の先頭と後尾の不連続により、 大きな折返し歪を生じ、正しい結果が得られなくなります。
  3. 不連続点の値として下記の中心値を使えば逆変換が成り立つ。
      f(t) = (f(t+) + f(t-)) / 2                   (3.1.3)
    
F(f) が周波数 fc で帯域制限されていれば、 標本化間隔を 1/(2*fc) に選ぶことで、 離散フーリエ変換が正しい近似になるわけですが、 この fc はナイキスト周波数(Nyquist sampling rate) と呼ばれます。

また、現実の非周期関数に離散フーリエ変換を適用する場合は、 後述のように、 元波形に乗算することで強制的に周期関数に変換して、 折り返し歪を減らすための窓関数(Window)がいろいろ 開発されています。

3.2. FFT と離散フーリエ変換

離散フーリエ変換では N^2 回の複素数乗算と N*(N-1) 回の複素数加算という多量の演算が必要になりますが、次式の計算で、 指数関数の周期性を活かした FFT(Fast Fourier Transform) 演算を使うと N^2 回の演算を N*log2(N) 回に減らすことができて、 計算時間は劇的に短縮されます。

         N-1
  F(n) = Σ(f(k)*exp(-j*2*π*n*k/N)  (n = 0, 1, .., N - 1)    (3.2.1)
         k=0
ただし、 FFT 演算の結果は離散フーリエ変換そのものではなく、下記の注意が必要です。
  1. データ数 N は 2 の冪乗(1, 2, 4, 8, 16, ..)に制限される。
  2. FFT 演算では、変換データの後半分に、 周波数軸負側のスペクトルを +2*(N-1) 移動した値が入るので、 実部は t = N/2 に対して線対称、 虚部は (N/2, 0) に対して点対称になる。(N/2 に於ける虚部の値は 0)
      実部 : 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)
  3. 時間軸のサンプリング間隔を T として、FFT 演算結果の実部を T 倍し、虚部を -T 倍すると離散フーリエ変換になる。
  4. 周波数軸のサンプリング間隔を F として、入力データの虚部の符号を反転してから FFT 演算を行い、出力データの実部を F 倍、虚部を -F 倍すると離散フーリエ変換にな る。(注5)

データ数が 2 の冪乗(1, 2, 4, 8, 16, ..)でなければばらないという制限については、元データの一部を削るか データの最後に 0 を追加して 2 の冪乗に合わせるのが普通ですが、 後述のステップ状関数補正を使う場合は最後の値を延長します。

また、データ数の制限を完全に回避する方法として、 chirp z-transform algorithm(Bluestein's FFT algorithm)(注6) があって、演算量は数倍に増えますが、用途の広い手法です。

3.3. 窓関数

非周期時間関数を離散フーリエ変換すると、 波形の終端が始端につながる周期関数と見られるため、 両端の不連続性に起因する大きな折返し歪を生じ、 正確なスペクトルが得られませんので、 元波形の両端が一致するような関数を掛けて 周期関数にしてしまう手法が考案され、 この種の関数を窓関数(Window)と呼んでいます。 なお、窓関数を掛けない状態は定数値 1 を取る窓関数を掛けたのと同じですから、 これを箱型窓関数(box-car window)と呼びます。

いずれにしても、波形の両側の振幅を除々に小さくして行って両端を繋ぐ仕組みですか ら、情報の一部が失われますし波形の面積も減ってパワーも減少しますが、折返し歪の 影響で識別できなかったスペクトルが分離できるようになるのが利点で、折返し歪の減 少を重視したもの、周波数分解能を重視したもの、パワー減少率を少なくしたものなど、 いろいろな種類があります。

始端と終端の不連続をなくす目的で、もっとも良く使われるのは、 下記のHanning窓で、

  x(t) = 0.5*(1 + cos(2*π*t/T))
  ここに、
	t = 時間
	T = 打切り幅
窓関数を使わない場合に比べて、 不連続波形の小さなスペクトルでも検出できるようになりますが、 パワーは 0.375 倍 (-4.26 dB) に低下します。

この他、大きな周波数成分近くの小さな周波数成分の分離に適した Hamming窓、

  x(t) = 0.54 + 0.46*cos(2*π*t/T)
	(パワー減少率 4.0 dB)
周波数分解能の良いBlackman-Harris窓、
  x(t) = 0.423 + 0.496*cos(2*π*t/T) + 0.0792*cos(4*π*t/T)
	(パワー減少率 5.14 dB)
線スペクトルの各周波数成分の比を正確に求めるのに適した flat-top
  x(t) = (0.54 + 0.46*cos(2*π*t/T))*sin(4π*t/T)/(4*π*t/T)
	(パワー減少率 7.00 dB)
など、いろいろあります。(注7)

4. ステップ状関数の離散フーリエ変換

現実の波形には、f(t) の値が 0 から始まって 0 以外で終わる ステップ状(step-like)の非周期波形がたくさんありま すが、

  f(t) = 0   (t <= 0)
       = f(T)  (T < t)
すなわち、ある時間 T 以降は変化しない場合は、時刻 0 から T までの波形を離散フ ーリエ変換した結果に、時刻 T から立ち上がるステップ波形のスペクトル
  f(T)/(1 - exp(-j*k*2*π/N))  (k = 1, 2, ..., N - 1)   (4.1)
を加算するだけで、正しい結果が得られることが わかっていて(注8)、 システムのステップ応答を調べる TDR (Time Domain Reflectmetry) 波形などの 解析で威力を発揮します。

波形が 0 以外から始まるときは、

  (f(T) - f(0))/(1 - exp(-j*k*2*π/N))  (k = 1, 2, ..., N - 1)   (4.2)
を加算します。

(4.1) を FFT 演算の結果に加算するとくは、FFT 演算の虚部の符号が反転している ため、虚部は減算しなければなりません。また、FFT でなく、DFT 演算の結果を補正 する場合は、下記のようになります。

  f(T)/(1 - exp(-j*2*π*f*dT))  (k = 1, 2, ..., N - 1)   (4.3)
  ここに、
	f = 周波数 (Hz)
	dT = 時間軸のサンプリング間隔 (s)

5. 離散フーリエ逆変換とラプラス逆変換

ステップ状波形など周期性のない時間関数のスペクトルを(周期性を前提とする)離散 フーリエ逆変換して時間関数に戻すと、 大きく歪んだ周期関数になってしまいますが、 (周期性を前提としない)数値ラプラス逆変換を使うと大幅に改善されます。

ラプラス変換とラプラス逆変換(注9)の定義は紛れがなくて、 それぞれ、

          ∞
  Fl(s) = ∫f(t)*exp(-s*t)*dt                  (5.1)
        -∞

                      g+j*∞
  f(t) = (1/(2*π*j))*∫Fl(s)*exp(s)*ds        (5.2)
                     g-j*∞
ですが、
  exp(-a*t)*f(t)*U(t) = ∫Fl(a + j*2*π*f)*exp(j*2*π*f)*df      (5.3)
が証明できますので(注10)、 ラプラス逆変換を下記のように計算することができます。
                   ∞
  f(t) = exp(a*t)*∫Fl(a+j*2*π*f)*df           (5.4)
                 -∞

(5.3) 式を見ると、f(t) に exp(-a*t) という強い減衰がかかっているため、周期関数 に近付いて、 離散フーリエ逆変換できない関数でもうまく扱えるの と、微分方程式の初期条件を自然に入れられるのが離散ラプラス逆変換の利点です。 ただし、a を大きくして収束を早めると後で乗算するとき発散し、 かといって a を小さくすると収束しなくなりますから、 適度な値を選ぶ必要があって、 5 <= a*T <= 7 (T = 打ち切り幅) ぐらいが適当という のが経験則です。

このメカニズムは、非周期関数に exp(-a*t) を掛けてフーリエ変換し、 その結果をフーリエ逆変換してから exp(a*t) 倍する計算をしてみると、 よくわかると思います。

この方式によるラプラス逆変換の手順は、次のようになります。

  1. ラプラス変換域の関数を F(s) として、 下記のデータ列を作成する。
      F(k) = F(a + j*2*π*k/T)/T  (k = 0, 1, 2, .. N/2)
      ここに、
    	a = 実定数 (5 - 7 程度)
    	T = 解が必要な時間幅
    	N = サンプリング数 (2 の冪乗)
    
  2. 上記のデータに離散フーリエ逆変換を実行する
  3. 得られた時間軸データを exp(a*k) 倍する (k = 0, 1, .. N)

離散フーリエ逆変換のプログラムで正規化(1/N 倍)していない場合は、 F(k) を N 倍しておかなければなりません。

例えば、時間 t0 から立ち上げるステップ関数

  F(s) = exp(-t0*s)/s
をラプラス逆変換して、いくつかの a の値を試してみると良いと思います。 このような帯域の広い関数では、折返し歪も exp(a*t) 倍されるため。 時間軸波形の後半で大きな振動が発生しますが、 Hanning などの窓関数を併用すると良い結果が得られます。 また、N を 2 倍にして、後半を捨てるといったやりかたもあります。

なお、この手法は回路解析でも有効で、電力分野で使われている EMTP (Electro Magnetic Transients Program) や ATP (Alternative Transients Program) などと同等の結果をより効率的に得ることができます。

この手法を使えば、下記の手順で、 非周期関数を離散フーリエ逆変換することができます。

  1. 周波数領域の関数を F(f) として、 下記のデータ列を作成する。
      F(k) = F(k/T - j*a/(2*π))  (k = 0, 1, 2, .. N/2)
      ここに、
    	a = 実定数 (5*T - 7*T 程度)
    	T = 解が必要な時間幅
    	N = (時間軸の)サンプリング数 (2 の冪乗)
    
  2. 上記のデータに離散フーリエ逆変換を実行する
  3. 得られた時間軸データを exp(a*k/N) 倍する (k = 0, 1, .. N-1) (時刻は k * T / N)
時刻は k * T / N になります。 実体は離散ラプラス逆変換と同じですから、 a や窓関数については離散ラプラス逆変換と同じです。

6. TDR/TDT 波形から S パラメータへの変換

TDR/TDT (Time Domain Refrectmeter/Time Domain Transmission) 測定は、 立上りの速いステップ波形を試料(DUT - Device Under Test) に印加して、その反射波形と透過波形を観測する単純な仕組みですが、 伝送路の特性インピーダンス、伝搬速度(位相特性)、減衰、不連続性、 インダクタンスやキャパシタンスなどの集中定数素子などが直観的に見られて、 非常に有効なツールです。 ただ、時間領域だけでなく、周波数領域で見たいことも多く、以下、 TDR/TDT の反射波と透過波の時間軸波形から試料の周波数特性を求めることを考えます。

6.1. ステップ応答とインパルス応答の関係

S パラメータなど周波数領域の特性は試料のインパルス応答ですし、 TDR オシロスコープは試料のステップ応答を見ますから、 ステップ応答からインパルス応答への変換が必要になりますが、 インパルス関数はステップ関数の時間微分ですから、 周波数領域で考えると簡単で、下記のようになります。(注10)


  インパルス応答のフーリエ変換 = j*2*π*f*{ステップ応答のフーリエ変換}
  ステップ応答のフーリエ変換 = {インパルス応答のフーリエ変換}/(j*2*π*f)
  ここに、
	j = sqrt(-1)
	f = 周波数 (Hz)

ここで、「線形システム」の章で扱った R-C 直列回路のステップ応答からシステム関数を求める例を見直してみてください。

6.2. ステップ応答から S パラメータへの変換

システムの応答を周波数領域で考える場合は、 S パラメータを使うことが多く、 測定にはベクトル・ネットワーク・アナライザ (VNA - Vector Network Anayzer) を使うのが普通ですが、高周波に於ける試料の基本的な特性

  S11(f) - 反射波 / 入射波    (6.2.1)
  S21(f) - 透過波 / 入射波    (6.2.2)
を電圧比で測定し、電力比のdB で表現するのが習慣です。 S11 は入力端から見た負荷インピーダンスを反映し、S21 は伝達関数になります。

TDR/TDT 波形から S11, S21 を求める場合は、 入射波(ステップ入力電圧)で規格化した反射波、 透過波をそれぞれフーリエ変換してから、 j*2*π 倍してインパルス応答に変換し、 (6.2.1), (6.2.2) を計算すればよいわけで、 下記のようになります(注12)。


  S11(f) = j*2*π*f*(反射波のフーリエ変換)    (6.2.3)
  S21(f) = j*2*π*f*(透過波のフーリエ変換)    (6.2.4)

例えば、VNA のシステムインピーダンスを R として、 Port1 に R-C 直列回路を接続して、S11 を求める場合を考えると、 反射電圧/入射電圧は

  ρ(t) = 1 - exp(-t/(R*C))
ですから、ρ(t) をフーリエ変換して、
  δ(f)/2 + 1/(j*2π*f*(1 + j*2*π*f*R*C))
j*2*π*f 倍すれば、
  S11(f) = 1/(1 + j*2*π*f*2*R*C)
一方、R-C 回路のインピーダンスは Z = R + 1/(j*2*π*f*R*C) ですから、 周波数領域の伝送理論で考えた反射係数は、
  S11(f) = (Z - R)/(Z + R) = 1/(1 + j*2*π*f*R*C)
で、時間領域から変換したものと一致します。

また、VNA のシステムインピーダンスを R として、 Port1 と Port2 の間にキャパシタンス C を挿入した回路を考えると、 透過電圧/入射電圧 τ(t) は C-R 直列回路に 1 V のステップ電圧を加えたときの R (Port2) の電圧ですから、

  τ(t) = exp(-t/(R*C))
τ(t) をフーリエ変換して、
  R*C/(1 + j*2*π*f*R*C)
j*2*π*f 倍して、
  S21(t) = R*C/(1 + j*2*π*f*R*C)
一方、周波数領域の交流理論で R (Port2) の電圧を求めると、
  S21(f) = R*C/(1 + j*2*π*f*R*C)
で、時間領域から変換したものと一致します。

(6.2.3), (6.2.4) の注意点は、 反射波も透過波も 0 から始まる波形にしておかなければならないことで、 TDR の場合は TDR ステップパルス波形の ρ == -1 から ρ == 0 になるまでの部分を 切り捨てるか 0 に書き換える必要があります。

6.3. システム誤差の校正

TDR/TDT システムに誤差がなければ、観測結果から反射波と透過波を分離して、 入射波との比を求め、 (6.2.3), (6.2.4) で変換すれば S11, S21 が求まりますが、 現実の測定器では、ソース側と負荷側の反射、伝送特性について、

  1. 方向性誤差 (Directivity - 反射波に入射波が混入)
  2. インピーダンス誤差 (Source Match/Load Match - システムのインピーダンス誤差)
  3. 周波数特性誤差 (Frequency Responce, Tracking - 周波数特性が平坦でない)
  4. ポート間相互(洩れ)結合 (Isolation - 透過波に入射波が混入)
といった誤差要因があって、精度の良いデータを得るには、 得られたデータの校正(calibration)が必要です。

VNAでは 2 ポート測定の場合、

  1. Short (短絡標準)
  2. Open (開放標準)
  3. Load (システムインピーダンス標準 - 通常 50 Ohm)
  4. Thru (ポート間直結標準)
の 4 つ標準を使って、 順方向と逆方向それぞれ、6 つの誤差要因を校正するのが普通です。 多少面倒な手続きですが、この手法こそが VNA 技術の基盤なのです。

TDR/TDT の場合も、得られた波形をフーリエ変換して、 周波数領域で VNA と同じ校正を行えば、より精度の高い S パラメータが得られますし、 周波数領域で校正した後、再度、フーリエ逆変換すれば、 校正された時間軸波形が得られます。

この校正についても、いろいろな手法が開発されていて(注13)、 以下、Short/Open/Load/Thru 標準を使った、いくつかの手順を説明しますが、 高い周波数では Short 標準の残留インダクタンスと遅延、 Open 標準の残留キャパシタンスと遅延、 全ての標準の損失と遅延が無視できなくなって(注14)、 標準に添付された正しい値を使わなければならなくなりますから、 多少繁雑になります。

6.3.1. 1 Port 補正

数 100 MHz 程度までの周波数で、反射係数 -1, +1, 0 という 純粋な Short/Open/Load 標準があれば、 下記の手順で、入力端子のみを持つ試料の S11 測定の 方向性誤差、インピーダンス誤差、周波数特性誤差を補正することができます。

実際の測定では、反射波の影響がなくなるまでの十分な時間幅が必要なことと、 反射波と透過波以外の波形を入れてはならないことに注意してください。 また、すべての測定でサンプリング間隔を同じにします。 波形の立ち上がり時刻はスペクトルの位相を変えるだけですから、 無理に一致させる必要はありません。

  1. TDR の Port1 に Open 標準を接続して、反射波形 v_open(t) を測定する
  2. TDR の Port1 に Short 標準を接続して、反射波形 v_short(t) を測定する
  3. TDR の Port1 に Load 標準を接続して、反射波形 v_load(t) を測定する
  4. TDR の Port1 に 試料(DUT) を接続して、反射波形 v_dut(t) を測定する
これらの時間軸波形をフーリエ変換して、下記のスペクトルを求めます。 電圧 0 から立ち上がる反射波と透過波だけを変換するために、 波形の立上り以前を 0 に置き換えるか削除して 0 から始まる反射波形にしなければなりませんし、 反射電圧を 0 から 1 の範囲に正規化するために、 ベースラインの変更やスケーリングが必要になることがあります。
  1. Vopen(f) = (v_open(t) のフーリエ変換)
  2. Vshort(f) = (v_short(t) のフーリエ変換)
  3. Vload(f) = (v_load(t) のフーリエ変換)
  4. Vdut(f) = (v_dut(t) のフーリエ変換)
このフーリエ変換ではステップ状波形に関する (4.1) の補正を使います。 データ数が 2^n に満たないときは、 最後の波高値を延長します。

試料(DUT)の正しい S11 を下記のように計算します。 すべて複素数です。

  S11 = (Vdut - D)/(S*(Vdut - D) + R)
  ここに、
	D = Vload
	S = (Vopen + Vshort - 2*Vload)/(Vopen - Vshort)
	R = 2*(Vopen - Vload)*(Vload - Vshort)/(Vopen - Vshort)
測定システムに誤差がなくて、Vshort = -1, Vopen = +1, Vload = 0 の場合は S11 = Vdut になることを確認してください。

高い周波数で、Short/Open/Load 標準のインダクタンス、キャ パシタンス、遅延、損失を考慮する場合は、D, S, R を下記のように求めます。

  S = ((Vopen - Vload)*(Sload - Sshort) - (Vload - Vshort)*(Sopen - Sload))
	/ ((Vopen*Sopen - Vloas*Sload)*(Sload - Sshort)
		- (Vload*Sload - Vshort*Sshort)*(Sopen - Sload))
  R1 = ((Vopen - Vload)*(Vload*Sload - Vshort*Sshort)
		- (Vload - Vshort)*(Vopen*Sopen - Vload*Sload))
	/ ((Sopen - Sload)*(Vload*Sload - Vshort*Sshort)
		- (Sload - Sshort)*(Vopen*Sopen - Vload*Sload))
  D = Vload - Sload*R1 - Vload*Sload*S
  R = R1 + D*S
  ここに、
	Sshort = Short 標準の反射係数
	Sopen = Open 標準の反射係数
	Sload = Load 標準の反射係数
Sshort, Sopen, Sload は標準インピーダンスの終端インダクタンス、 終端キャパシタンス、遅延、損失を考慮した値ですから、複素数になります。

6.3.2. 2 Port 補正

入力端子と出力端子を持つ 2 Port 試料を TDR/TDT の Port1 と Port2 の間に挿入し、S11, S12 を測定する場合は、下記の手順で、すべての 誤差を補正します。

  1. 前記の手順で、D, S, R を求める
  2. Port1 と Port2 の両方に Load 標準を接続し、v_isolation(t) を測定する
  3. Thru 標準で Port1 と Port2 を直結し、Port1 の波形 v11_thru(t) と Port2 の波形 v21_thru(t) を測定する
  4. 試料を Port1 と Port2 の間に挿入し、Port1 の波形 v11_dut(t) と Port2 の波形 v21_dut(t) を測定する
  5. 試料の Port1 と Port2 を反転し、Port1 の波形 v22_dut(t) と Port2 の波形 v12_dut(t) を測定する
これらの時間軸波形をフーリエ変換して、下記のスペクトルを求めます。
  1. V_isolation(f) = (v_isolation(t) のフーリエ変換)
  2. V11_thru(f) = (v11_thru(t) のフーリエ変換)
  3. V21_thru(f) = (v21_thru(t) のフーリエ変換)
  4. V11(f) = (v11_dut(t) のフーリエ変換)
  5. V12(f) = (v12_dut(t) のフーリエ変換)
  6. V22(f) = (v22_dut(t) のフーリエ変換)
  7. V21(f) = (v21_dut(t) のフーリエ変換)

ここでも、 TDR/TDT 波形を 0 から 1 の範囲に正規化するためのデータのスケーリングやベースラ イン変更を忘れないでください。例えば、TDT 波形の電圧は、-250 mV から 0 V の範囲 ですから、250 mV を加算して、4倍にする変換が必要です。

試料(DUT)の正しい S11, S12 を下記のように計算します。

  S11 = (1 + (V22 - D)*S/R)*(V11 - D)/(R*P)
		- (V21 - X)/T*(V12 - X)/T*L/P
  S12 = (1 + (V11 - D)*(S - L)/R)*(V12 - X)/T/P
  ここに、
	X = V21_isolation
	L = (V11_thru - D)/(R + S*(V11_thru - D))
	T = V21_thru * (1 - S*L)
	P = (1 + (V11 - D)*S/R)*(1 + (V22 - D)*S/R)
		- (V21 - X)/T*(V12 - X)/T*L^2
通常 X は極めて小さいので 0 と看做すことができます。 P は 1 に近い値になります。

V22, X12, V21 は試料に方向性がある場合、試料の Port1 と Port2 を反転させて測定 するか、測定器にその機能があれば、Port2 側からステップ波形を加えて測定します。

試料に方向性がなければ、 V11 = V22, V12 = V21 として計算し、 S21 = S12 とすることができます。

6.3.3 校正の簡略化

通常の測定器は高いシステム・インピーダンス精度を持っていますから、 インピーダンス誤差を無視すれば、 下記のように簡略化できます。

  S11 = (V11 - Vload)/(Vshort - Vload)
  S21 = (V21 - V_isolation)/(V21_thru - V_isolation)
前者は Short-Load 補正、 後者は V_isolation を無視すれば Thru 補正になっていますが、 周波数が低ければ、この程度の校正で間に合いますし、 通常は V_isolaton = 0 と看做すこともできます。

6.3.4 帯域と分解能

(3.1.1) から

  Bw = 1/(2*dT)
  T = 1/dF
  ここに、
	Bw = 帯域幅 (最高周波数) (Hz)
	dT = サンプリング間隔 (s)
	T = 時間幅 (s)
	dF = 1/(N*dT) (周波数領域のサンプリング間隔)
の関係があって、 一般の TDR/TDT では、サンプリング間隔 10 ps、データ点数 512、 時間幅 5 ns 程度ですから、 帯域幅 50 GHz、周波数間隔 200 MHz 程度のデータが得られるはずですが、 オシロスコープの帯域が 20 GHz 程度ですから、 これが上限になります。

逆に、VNA の周波数データをフーリエ逆変換で時間軸に変換する場合も 同じような限界があって、例えば、20 GHz 帯域で 512 点の測定ができたとしても、 周波数間隔 40 MHzですから、 時間間隔 500 ps、時間幅 25 ns の時間軸データしか得られません。

TDR/TDT の波形から S パラメータを求める場合、 4 - 6 GHz 程度なら校正なしとか、 簡略化した校正で間に合いますし、 Short-Open-Load-Thru (SOLT) 校正をすれば 8 - 12 GHz まで使えるのが普通です。

ケーブルのような長い伝送路の測定では、 反射波の影響がなくなるのに長い時間がかかりますから、 サンプル数を 4095 点程度に増やしたとしても、 サンプリング間隔を長くせざるを得なくて、 帯域が狭くなります。

なお、VNA は狭帯域で長い時間をかけた測定を行えば、 広いダイナミックレンジが確保できますし、 測定対象も TDR/TDT より広いのですが、 TDR/TDT は直観的で測定も簡単というのが利点で、 TDR/TDT を使った S パラメータの測定にも捨てがたいところがあります。

7 注

7.1. フーリエ変換の定義

フーリエ変換には、ここで引用した文献でも、Papoulis では下記の定義、 Brigham はこの解説の定義と、係数の違いによる、いろいろな可能性があって、 下記の式の定義だと、システム関数に 2*π が含まれないため、数式が簡潔になる半面、 実務で使われる周波数を考えるとき ω から f への換算が必要になります。 この解説でも、この下記の定義のほうが数式が簡単になる部分が多いのですが、 実務で周波数がほしいことが多く、前記の定義にしてみました。

  Harry Bateman,- TABLES OF INTEGRAL TRANSFORMS
	(McGRAW-HILL) 1954
でも、これを使っていますが、このテキストを書き始めてから、 この解説でも下記の定義にすべきではなかったかという疑問が残りました。
          ∞
  F(ω) = ∫f(t)*exp(-j*ω*t)*dt
        -∞
       = R(ω) + j*X(ω)
       = A(ω)*exp(j*θ(ω))

                    ∞
  f(t) = (1/(2*π))*∫F(ω)*exp(j*ω*t)*dω
                   -∞
  ここに、
	t = 時間 (s)
	ω = 角周波数 (rad/s)
	   = 2*π*f
	f = 周波数
	π = 3.14159265..
	f(t) = 時間を変数とする関数
	F(ω) = f(t) のフーリエ変換
	R(ω) = H(ω) の実部
	X(ω) = H(ω) の虚部
	abs(F(ω)) = sqrt(R(ω)^2+X(ω)^2) = 振幅 (スペクトル)
	θ(ω) = atan(X(ω)/R(ω)) = 位相角 (rad)
	j = sqrt(-1)
この定義とこの解説の定義ではフーリエ変換の結果に下記の違いがあります。
  H(f) = F(ω) / (2 * π)

7.2 時間不変な線形システムはインパルス応答のみから決定できる

例えば、Papoulis の (5-10) 参照してください。 線形性と不変性と言う極めて一般的な条件だけで、この驚くべき結果が出ます。 不変性というのはシステムの決定要因が時間的に変化しないということで、 システムが微分方程式で決められる場合は、その係数が定数になるという意味です。

7.3 時間不変な線形システムの固有関数と固有値

これも驚くべき結論ですが、 例えば、Papoulis の (5-11) 参照してください。

  f(t) = exp(j*2*πt) なら、出力は k*f(t)
ということで、 この比例定数 k は複素数ですが、時間不変な線形システムの固有関数が指数関数、 比例定数が固有値になって、 これが交流理論による回路計算の基盤になります。

7.4 R-C 直列回路のステップ応答

R-C 直列回路にステップ関数の電圧を印加したとき、回路に流れる電流は、

  (1/C)*∫i(t)*dt + R*i(t) = e(t)
  ここに、
	C = キャパシタンス (F)
	R = 抵抗 (Ohm)
	e(t) = R-C 直列回路に印加する電圧 (V)
	i(t) = R-C 直列回路に流れる電流 (A)
	i(0) = 0  (初期条件)
の簡単な線形微分方程式で、電気工学の「過渡現象論」だと、 素直に微分方程式を解くか、ラプラス変換で解くのが普通ですが、 i(t) のフーリエ変換を F(f) として、両辺をフーリエ変換すれば、
  F(f)*(1/(j*2*π*f*C) + R) = δ(f)/2 + 1/(j*2*π*f)
F(t) を求めると、
  F(t) = (1/R)*(1/(1/(C*R) + j*2*π*f))
フーリエ逆変換して、
  f(t) = (1/R)*exp(-t/(R*C))

7.5. FFT による離散フーリエ逆変換

複素数 z の共役複素数を <z>、つまり、<x+j*y> = x-j*y とすれば、 下記が成り立つので、逆変換にも同じプログラムが使えます。

  N-1                             N-1
  Σ(H(k*T))*exp(j*2*π*n*k/N) = <Σ<H(k*T)>*exp(-j*2*π*n*k/N)>
  k=0                             k=0

7.6. Bluestein's FFT algorithm

離散フーリエ変換

        N-1
  X(k) = Σ(x(n)*exp(-2*π*i*n*k)   (k = 0, 1, ..., N-1)
        n=0
  n*k = -(k - n)^2/2 + n^2/2 + k^2/2
を代入して
                          N-1
  X(k) = exp(-π*i*k^2/N)*Σ(x(n)*exp(-π*i*n^2/N)*exp(π*i*(k-n)^2)   (k = 0, 1, ..., N-1)
                          n=0
と変形すると、
              N-1
  X(k) = c(k)*Σ(a(n)*b(k-n))   (k = 0, 1, ..., N-1)
              N=0
  ここに、
	a(n) = x(n)*exp(-π*i*n^2/N)
	b(n) = exp(π*i^n^2/N)
	c(k) = exp(-π*i*k^2/N)
つまり、a(n) と b(n) の畳込積分に定数ベクトル c(k) を乗算するという構造で、b(n) と c(n) は共用できますから、b(n) の計算は1度だけ、 a(n) の計算と a(n)*b(k-n) の逆変換を毎回行えば済みます。

そして、これらの計算は FFT と同じ形をしていて、密度が半分ですから、 a(n), b(n) の後半に 0 を追加して 2 倍のデータ数に拡張すれば 任意のデータ数で FFT が使えるというアイデアです。

FFT 演算が 2 倍増えてデータ量も 2 倍になりますが、 FFT の高速性を活かしながらデータ数の制限からは開放されることになります。

  Leo I.Bluestein,- "A linear filtering approach to the computaion of the
	descrete Fourier transform", Northeast Electronics Research and Engineering
	Meeting Record 10, 218-219 (1968)
がこのアイデアですが、後に z 変換として一般化され、 chirp z-transform algorithmとしてよく知られるようになりました。

7.7. 窓関数

例えば、

  FREDERIC J.HARRIS,- On the Use of Windows for Harmonic Analysis with the
	Discrete Fourier Transform
	(PROCEEDINGS OF THE IEEE, Vol.66, No.1,January 1978)
にたくさんの評価例があります。

7.8 ステップ状関数の離散フーリエ変換

  A.M.Nicolson,- FORMING THE FAST FOURIER TRANSFORM OF A STEP RESPONSE
	IN THE TIME-DOMAIN METROGY
	(ELECTRONICS LETTERS 12th July 1973 Vol.9 No.14, pp 34-36)
がわかりやすいと思います。 これは、原波形から原点と波形が変化しなくなる点の波高値を通る傾斜関数を原波形か ら引いて、波形が変化しなくなった点で打ち切った周期関数の離散フーリエ変換値が、 周波数 0 以外で、原波形の打ち切り点までの離散フーリエ変換値に等しい(!)ことを利 用して、ステップ状関数の離散フーリエ変換を 0 で始まり 0 で終わる周期関数の離散 フーリエ変換に転化し、その結果に解析的に計算できる打ち切り点から立ち上がるステ ップ関数のスペクトルを加算することで、無限大まで伸びるステップ状関数のスペクト ルを求めるというアイデアです。

この他にも、1 標本区間ずらせた波形を引くとか、1 標本区間毎の波形の差を作って、 非変化点以降を 0 にするアイデアがあって、最終的に、これら3つのアイデアが、そ の発想の見掛けの違いにかかわらず、原理的に同じであるという解析結果が下記の論文 です。

  Jurg Waldmeyer,- Fast Fourier Transform for Step-Like Functions:
	The Synthesis of Three Apparently Different Methods
	(IEEE TRANSACTIONS ON INSTRUMENTATION AND MESUREMENT,
	Vol.IM-29,No.1,MARCH 1980, pp 36-39)

また、ステップ状関数の高周波部分の折返し誤差を減らすアイデアもいろいろあって、 例えば、

  G.D.Cormack and J.N.McMullin,- Frequency De-Aliased FFT Analysis of Step-Like
	Functions
	(IEEE Transactions on Instrumentation and Measurement, Vol.40, No.4, Aug 1991
では、FFT の結果に
  (sin(π*k)/(π*N))^2
を乗算するだけで、パワー計算の結果を少し改善してみせています。インタネットで公 開されていますので、ご覧になってみてください。例示された波形が面白いです。

7.9 ラプラス変換

s = a + j*ω と変換すると、ラプラス変換を exp(-a*t)*f(t) のフーリエ変換とみなす こともできて、FFT 演算のラプラス変換への応用に道が開けます。

ヘビサイド(Heaviside)の演算子法という天才のひらめきの帰結とも言える、ラプラス変 換にはたくさんの解説書があって、多分、最もよく知られたのは

  G.Doetsch,- Theorie und Anwendung der Laplace-Transformation
	(Dover, 1943)
だと思いますが、 学生時代に読んだ「実用ラプラス変換」(森北出版)という実務者向けのテキストが 明解かつ実用的でよかったです。

なお、この変換そのものは、ほぼ1世紀前のオイラーの発明で、Laplace の著作で 広く知られるようになったようです。 しばらく前に、オイラーの著作のごく一部が日本語で読めるようになりましたが、 今読んでも大きな感銘を受けます。

Heaviside の伝記の和訳がないのが残念。

7.10 ラプラス逆変換とフーリエ逆変換の関係

例えば、Papoulis (9-17) を見てください。 また、(5.2) のラプラス逆変換を

  s = a + j*2*π*f
で変数変換してフーリエ逆変換と比較しても良いと思います。

7.10 - ステップ応答とインパルス応答の関係

周波数領域では微分は j*ω 倍、積分は 1/(j*ω) 倍と代数化されますから、 ステップ関数の微分がインパルス関数になることに気付けば簡単ですが、 下記のように直接確認することもできます。

システム関数を H(f) = R(f) + j*X(f) とすれば、 ステップ応答のフーリエ変換はステップ関数のフーリエ変換を乗じた

  [δ(f)/2 + 1/(j*2*π*f)]*[R(f) + j*X(f)]
        = R(0)*δ(f)/2 + X(f)/(2*π*f) - j*R(f)/(2*π*f)
ですから、j*2*π*f 倍すれば、f*δ(f) = 0 に注意して、
  j*2*π*f*[δ(f)/2 + 1 / (j*2*π*f)] * [R(f) + j*X(f)]
        = R(f) + j*X(f)
        = H(f)
が得られます。f*δ(f) の値は、
   ∞
  ∫f(x)*δ(x)*dx = f(0)
 -∞
で f(x) = x とすれば、
   ∞
  ∫x*δ(x)*dx = 0
 -∞
から
  x*δ(x) = 0
で 0 になります。

7.11 - dB 表現

dB 表示を使う場合は

  -20*log10(S11(f)) = -20*log(S11(f))/log(10)
  -20*log10(S21(f)) = -20*log(S21(f))/log(10)
で変換してください。

7.12 - S11, S21 の計算

入射波がステップ関数なら、0 < f で、

  j*2*π*f = 1 / (入射波のフーリエ変換)
ですから、例えば、
 S11(f) = (反射波のフーリエ変換) / (入射波のフーリエ変換)
        = j*2*π*f*(反射波のフーリエ変換)
 S21(f) = (透過波のフーリエ変換) / (入射波のフーリエ変換)
        = j*2*π*f*(透過波のフーリエ変換)
です。

7.13 - 周波数領域に於けるシステム誤差の補正

いろいろありますが、VNA の場合は、例えば、

  J.Williams,- Accracy enhancement fundamentals for vector netwaork analyzers
	(Microwave Journal, Martch 1989, pp 99-114)

TDR/TDT の場合は、例えば、

  Tom Dhaene, Luc Martens, and Daniel De Zutter,- Calibtration and Normalization
	of Time Domain Network Analyzer Measurements
	(IEEE Trtansactions on Microwave Theory and Techniques, Vol.42, April 1994)

7.14 - Short-Open-Load-Thru 標準

Short や Open は特性インピーダンスの不連続面があって、 それが Short ならインダクタンス、Open ならキャパシタンスに見えます。 実際の値は、それぞれ、2e-12 H, 5e-15 F 程度ですが、 周波数の多項式で近似した計算結果が、 標準インピーダンスの製品に添付してあります。 また、標準インピーダンス素子から接続面までの伝送路(1 cm 程度)に於ける遅延と 減衰がありますので、すべての標準について、 インピーダンスの値とともに、 ケーブル延長で評価した損失と遅延の値が定義されていて、 標準インピーダンスは

  (定義値の)終端インピーダンスで終端された(定義値の)伝送線路の入力インピーダンス
として計算します。 Thru は直結なので、損失と遅延のない線路で評価できるのが普通です。

例えば、手元の 3.5 mm Calibration Kit だと、下記のようになっています。

ValueShortOpenLoadThru
LLs000
C0Co00
R00500
Offset Delay (ps)31.80829.24300
Offset Z0 (Ohm)50505050
Offset Loss (Ohm/s)2.36G2.2G2.2G2.2G
Ls = 2.0765e-12 - 1.0854e-22*f + 2.1705e-33*f^2 + 1e-44*f^3 (H)
Co = 4.943e-14 -3.1013e-25*f + 2.317e-35*f^2 - 1.6e-46*f^3 (F)
f = Frequency (Hz)

Offset Delay は位相遅れを速度係数 1 の線路の遅延時間で評価した値、 Offset Loss は遅延時間あたりの損失を 1 GHz に於ける抵抗で評価したもので、 物理的には導体の表皮効果になります。 等価線路の長さは 1 cm 程度ですが、1 GHz 程度の周波数では無視できなくなって、 インダクタンスとキャパシタンスが入れ替わってしまいます。

上記の定数から伝送線路定数への変換は、次のようになります。

  α*l = Offset_loss * Offset_delay / 2 / Z0 * sqrt(f/1e9)
  β*l = 2*π*f*Offset_delay
  ここに、
	α = 伝送線路の減衰定数 (neper/m)
	β = 伝送線路の位相定数 (rad/m)
	l = 伝送線路の長さ (m)
ほとんどの場合、α*l は無視できます。

平林 浩一, 2009-01-04