AD変換器の限界 その２

　高級品のAD変換器とワードクロック発振器を購入しました。

　AD変換器（オーディオインタフェース）は

1. Steinberg社製UR-824  http://japan.steinberg.net/jp/products/hardware/ur_series/lineup/ur824.html（７万円台）　
2. ワードクロック発振器は、長野県のメーカー,城下工業株式会社製のSWD-CL10. 　http://www.shiroshita.com/sw/swd-cl10.html　（６万円）

　AD変換器のマスタークロックは「ワードクロックジェネレータ」と呼ばれる機材で供給することができます。同社の製品にはTCXO級の水晶が搭載されており（オプションでOCXOにアップ可能）スペック上は0.28ppm となっています。

　早速試験してみましたが、AD変換器の内臓クロックでは10ppm程度外れているもの安定度は申し分なし。ワードクロック内臓発振器(TCXO)で動作させて試験した結果は、測定値で0.1ppm未満で今まで使用していた安物オーディオインタフェースとは段違いです。ただしこの水準にくると測定開始してからドリフトが目立ちます。

　最後にルビジウム発振器を接続して試験ですが、機材は通販で届いたばっかり。まだ冷え切っています。作業場のエアコン入れっぱなしにして一週間程度走らせっぱなし。すこしずつ安定してきて、チャンピオンデータで441.0 Hz の信号の測定値は, 441.000,000,000 ±50 (nHz)程度。周波数不確かさ換算で, 1.1 ×10^(-10) となりました。おそらくRb発振器の安定度限界か、FGとクロックジェネレータに使用されているDDS(Direct Digital Synthesizer)の限界と思われます。50Hzでも同様な試験をした結果も、ほぼ同じ水準。

　今度は、AD変換器のCH1,CH2に「同じ信号」を入力し、CH1,CH2間の位相差を測定。理想のAD変換器ならば、位相差=0が測定されるはずです。しかし結果はΔφ=50 μrad の位相差。 CH1,CH2のケーブルを逆にして同様に測定すると Δφ=60 μradの位相差

　これは何を意味するか？

　ケーブル長のわずかな差での遅延が影響しているのです。また定常的なAD変換器内での位相差は、おそらく電子部品のCH間での微妙な差が影響していると思われます。ちなみに信号レベルを大きくするとこの機器内での定常位相差は大きく変化します。この試験での10μradの差異がケーブル長由来でかつ光速の限界の影響だとすれば、441Hzの信号における1mの線路長の差は約10μradになります。（波長を計算してΔφ=L/λ ・2πと求められます) 　厳密に実験するには、インピーダンスによる遅延の影響をうけない「光伝送」で信号を送受信してみると面白いかもしれません。純粋に距離による位相差が観測されるはずです。。。

 世の中には100万分の1という精密なものがたくさんありますが、 このppmを超えると途端に様々な要因による誤差が表面化してきます。今回の測定値はppmを超えppbまで届きましたが、Rb発振器の不確かさ接近した水準であるため、これ以上の精度向上はむずかしいでしょう。

　ノイズの影響の軽減をねらって、バンドパスフィルタ(BPF)を前置していたのですが、実はBPFには測定揺らぎに影響する大きな落とし穴があったのです。次回以降この話題を取り上げたいと思います。

 

三角関数の変な公式？

　必要あって、次の計算式に取り組んでいました。

cos^2(θ) - cos^2(2θ) = sinθ・sin3θ　

cos^2(0) - cos^2(θ) =  sinθ・sinθ　（単純には 1-cos^2(θ)=sin^2(θ))

これを一般化して、

cos^2(nθ) - cos^2{(n+1)θ} = sin θ・sin (2n+1)θ　　を、とりあえず証明しました。

 

実は以前に説明した「リサージュ外積」 Lsをz領域で表したとき遅延段数kを1,2,3....としたときの一般形を知りたかったのです。 

展開と導出はそれほど難しくはなく、与式をまず「和差の積」に変換、次に三角関数の「和積変換」、sin cos の項に集めると「倍角の公式」を用いてさらに簡単になり目的の式ができあがりです。

　当社の周波数測定方法（計算方法）は、「リサージュ外積比(=r)」を解いて、目的の周波数を計算しているのですが、現在までにわかっている計算式は、1) r=Ls2/Ls1 とする「倍角法」、2)r=Ls3/Ls1 とする「三倍角法」　だけです。 

r = Ls5/ Ls3 のような形ではどうなるか・・・

結局、上の・・・"余弦自乗差・正弦倍角変換？　をつかうと、 r = sin(mΩ) / sin(nΩ)　という式で、与えたm,nの条件下で求めたい周波数Ωを解くことになるのです。　他に３つ以上のリサージュ外積を組み合わせた解法も（数学的には）可能ででしょうが、どんな性質や性能となるかは予想がついておりません。狙いとしては, Ls1,Ls2,Ls3 ....と適当なところまで使用した計算式。複数を組み合わせることでお得な何かがあればいいかと思っています。

　

ＡＤ変換器の限界 その１

 デジタル時代の盲点・・・AD変換器

　当社の方法をつかうと限られたサンプル数で非常に正確に「短時間周波数」を求めることができます。f = fs/2π*arcsin(0.5 * √(3- Ls3/Ls1) )　の方法を用いると合計５サンプルで計算できます。

　さて今回の話題ですが、以上のような数値計算の制約による計算誤差を工夫したとしてもある意味避けられない、アナログ信号をデジタル信号に変換する「ＡＤ変換器」の問題について説明したいと思います。

 　上の式にはサンプリング周波数fsが入っています。オーディオ用のＡＤ変換器（オーディオカード）のサンプリング周波数は普通44.1kHzが用いられます。

　では実機で動作しているサンプリング周波数がずれている場合、何がおこるか？

　ＡＤ変換器の内部では「水晶発振子」が使用され、サンプリング信号を生成していますが、実はこの水晶発振子がクセものなのです。水晶発振子の周波数は必ずしも完璧ではなく、数ppmから100ppm程度の誤差があります。パソコンクロック用など、周波数精度が問題ない場合は周波数精度はそれほど気にされません。

　水晶発振子には様々なクラスがあって、温度補償がついているTCXO、恒温槽はついて温度制御まで行っているOCXOなどのがあります。OCXOでは精度は0.02ppm程度となるものもあります。

　ここで「精度」と簡単に書いてしまいましたが、細かくは 1)公称値とどのくらいずれているか、2)動作中どのくらいの変動があるか、　の２点が重要です。

　さてここで問題発生。　試験に使用する「正しい周波数」をどう用意するか。ＡＤ変換器のサンプリング周波数は, 44.1 kHzからどのくらい外れているのか、どこまで信用して良いのか。

　実際に使用した「正しい周波数」は「ルビジウム発振器」の信号を外部同期入力とした信号発生器(Function Generator=FG)を使用しています。ルビジウム発振器の「精度・安定度」は"数十ppb"クラスとされています。

　実際に測定してみると・・・

　測定値が、常に -35～ -40ppm くらい低いのです。50.000,000 Hzを測定すると、49.998 ～49.999Hzの間となっています。エージングという慣らし運転を繰り返しＡＤＣを安定化すると、誤差は-39ppm～-41ppmに落ち着いてきました。搭載されているデバイスは普通のものなのでしょう。

　当社の方法を用いて高精度な測定のためには、正しい値のわかっている信号を用いて、以上説明したような「ＡＤ変換器の校正」が必要になります。

　資金があれば、「ワードクロック同期」のついたＡＤ変換器を使用、さらにワードクロック発生器もOCXOクラスの高級品を使用すれば精度（不確かさ）は向上します。

　前述の「安定度」ですが、測定結果が, 20μHz程度のドリフトと,数秒間での測定値の不確定さは約10μHzとなっています。20μHzの測定値のズレは0.4ppmに相当し、公称値が40ppmズレているとはいえ変動はそこそこであることがわかっています。

　かつて大学院在学中での「旧方式」による試験（※ 本ブログで書いてあるものとは若干異なります）でとれた「チャンピオンデータ」では、441.000,000 Hzの試験信号（当然 Rb発振器でロック）、ADCやクロックジェネレータも最高級品を使用し、1μHz級の周波数分解能、位相差測定では1μradの分解能が確認できました。　この検証については論文で公表しておりますので興味のあるかたはどうぞ。

文献：　西江, 赤木, 弦楽器F0 推定のための精密周波数測定方法, 電子情報通信学会論文誌 A J97-A(4), 332-342, 2014-04-01.  http://hdl.handle.net/10119/11937

ーーーー

　Rb発振器で同期した環境を用いると、分単位の測定値ではおそらく有効数字8桁は実現できていると思います。しかし「なぜ８桁なのか」？

　実は・・・AD変換器の分解能は16ビット。当社の方法は振幅方向の分解能も周波数分解能に関係しているのです。シミュレーションのみでしか検証はできていませんが、「リサージュ外積Ls」の分解能の限界と考えています。クロス積の計算では「桁落ち」が発生してしまうのです。

 

 

 

DFT/FFTの限界 ー不確定性ー

DFT/FFT の限界

DFT/FFTを用いた周波数測定の限界について説明したいと思います。

DFT/FFTのサンプル数とサンプリング周波数に応じて、「周波数ビン」とも呼ばれる目盛りが決定されます。50Hzの測定で仮に1Hzの分解能が欲しければ、サンプル数Nとサンプリング周波数fsの間には以下のような条件が想定されます。

・　N=4096, fs=4096 (Hz)　

　こうすると、測定に要する時間は１秒となります。では測定時間分解能を0.1秒にしたい場合どうすればよいか。サンプリング周波数を高くする方法をまず考えてみます。

・　N=4096, fs =40960 (=40.96kHz) 

すると周波数目盛りである周波数ビン＝10Hzとなってしまいます。一見周波数分解能が低下するように見えるのですが、得られたスペクトルF(ω)の実部・虚部をベクトルに見立てると回転が起こっており、この回転率をみることで、周波数ビンの分解能を超えた計測が可能です。

　仮に、１測定あたり、スペクトルの角度が36度変化しているなら、本条件の場合、1/10秒で36度の回転なので、１秒あたり360度、すなわち基準となる周波数ビンから1Hzの差があることがわかります。測定限界は180度です。これを超えるとベクトルの回転が右回りなのか左回りなのか判別できません。180度/0.1秒の変化は5Hzに相当します。つまり、隣接する周波数ビンの間は、前述のようなベクトルの回転をふくめればなんとか計算できることになります。

　逆にサンプル数Nを1/10にしても同じことで（FFTでは２の倍数の制約があるので現実的ではありませんが）、周波数ビンの間隔が10倍になる代わりに測定時間分解能は1/10になります。

　まとめると、「周波数ビン」という「目盛り」を一定にしたいなら、『測定時間（時間分解能）とサンプリング周波数は反比例の関係』にあり、サンプリング周波数を一定にするなら、周波数ビン（＝目盛り）と測定時間（時間分解能）はこれも反比例の関係。

　DFT/FFTでいう周波数ビンの間隔＝周波数分解能を高くすると、測定時間が長くなるジレンマが存在するわけです。

　一般的には「フーリエ変換の不確定性」と呼ばれていますが、不確定性の厳密な定義はちょっとややこしいです。時系列波形s(t), スペクトルのS(ω)　の「１次モーメント」（＝平均持続時間、平均周波数）と「２次モーメント」（＝平均時間窓、周波数分解能）で定義されます。　参考： L.Cohen／吉川、佐藤訳　時間-周波数解析　（朝倉書店）　

　以上サンプリング周波数fsとサンプル数Nの関係で説明したように、DFT/FFTのような離散系の例ではコーエンのような面倒な積分計算なしに、なにかしら周波数分解能と時間分解能のあいだにトレードオフがあることは、感覚的に理解できるでしょう。

　もうひとつ、DFT/FFTを使用する際には考慮しないといけない事項があります。DFT/FFTでは最低１サイクルの信号がないといけません。周波数ビン=1の位置にある周波数が周波数の下限となります。

　このことをふまえ電源の系統周波数(50,60Hz)を精密測定することを考えてみましょう。fs=44100 Hzという一般的なオーディオ用ADCを使用し、測定時間を概ね0.1秒にするなら、 N=4096が候補となります。周波数ビンは10.77Hzです。さらに一歩つっこんで、1mHzの周波数分解能を求めるなら、複素スペクトルの回転角度はおおまかに 1/10000回転=0.036度となります。

　はたして、得られたスペクトルからこの角度速度の計算は現実的な水準か。DFTを使用する1mHz on 50Hzの測定では技術的な課題となります。

　周波数スペクトル（複素数）の「回転」をどこまで精度良く求められるかが、周波数測定精度を決定してしまいます。次回は「複素ベクトルの回転の計算方法」について書きたいと思います。

一般的に知られる方法：スペクトルの実部・虚部から φ=arctan ( Im / Re ) と角度を求めておき、その角度の差分をとる以外にエレガントな方法があります。