TNJ-008：正弦波をA/D変換し窓関数なしに打ち切ってFFTしてみると

 2014年1月20日 公開

はじめに

こんなご質問をいただきました。

「A/Dコンバータのデータシートに記載されているFFTスペクトルについてなのですが…」

「A/D変換サンプリング周波数fS [Hz]（サンプリング周期TS [s] = 1 / fS）とFFTポイント数Nから決まるFFT時間長TFFT [s]（TFFT = N / fS）と、入力信号（被測定信号）周波数fIN [Hz]の周期TIN [s]の整数倍の時間（TMES [s] = PTIN, Pは任意の整数」がぴったり合わないようなケース（TFFT ≠ PTIN）があると思います」

このご質問はサンプリングが正弦波周期の途中で「窓関数なしに」打ち切られると（正確にはサンプリングされた波形のオシリとアタマが連続になっていない）おかしくならないか？ということを意図しているのでした。ご質問は続きます。

「しかしA/Dコンバータのデータシートでは、窓関数を用いていないようなスペクトルになっています。たとえば（図1に示す）AD9626のデータシートFig. 11の条件なのですが…」

「データシートはほとんどこのケースと思います。これはどう考えればよいのでしょうか？」

結果的にこのご質問の答えとしては、この技術ノートの検討のように「サンプリングが正弦波周期の途中で打ち切られると、スペクトルはおかしくなる」ということなのでした。

図1. A/Dコンバータのデータシートに記載されているFFTスペクトルの例（AD9626）

FFT時間長と被測定信号の波数（サイクル数）がぴったり合わないときはどうなるの？

ということでご引用いただいたAD9626は…

AD9626: 12ビットA/Dコンバータ、1.8V、170/210/250 MSPS

またそのFig. 11の条件とは、

サンプリング速度 = 170MSPS

信号周波数fIN = 10.3MHz（シングル・トーン）

FFTポイント数 = 64k（65536ポイント。216）

この設定（FFT時間長 TFFT = 385.5059us）であれば、FFTを行う期間内に10.3MHz信号が3970.7…波数（サイクル数）相当になります。これを3971に波数（サイクル数）として丸め、64×1024（= 65536）×10.3MHz / 3971=169.9876102MHzと計算し、得られたものをA/D変換サンプリング周波数fSとすればよい…、ということなのか？というご質問でありました。

 

普通は窓関数を掛け合わせるのだが

普通はこのような（信号周期を考えない）場合は、信号自体に窓関数を掛けて、それにより暴れを無くします。これはよく聞く話ではないでしょうか。私は「これはギブス現象に関連するのではないか…」と思い、「おっしゃるように窓関数を用いていなければギブス現象により、周波数スペクトルが暴れることがありますね」と返しました。その続きとしても以下をコメントしました。とはいえこれは「間違い」でありました…（そのため、以下の節は読まずにスキップしてください）。

 

嗚呼かんちがい…のコメント

「このプロットにはいくつか理由が考えられるかと思います。データがナイキスト周波数まで出ているので、おっしゃるようなギブス現象として見えるのではないかと感じられますね。

1. ギブス現象は信号のスペクトルとサンプリング長に相当する矩形フィルタとの周波数軸の畳み込みf(ω) ⊗ sinc(ω)となるはずです（⊗は畳み込み演算を表す。畳み込みは追って示す）。ここでは信号がシングル・トーンとなっていますので、それが目立たないのではと思われる点
2. それではノイズ・フロアがうねる（畝る）はずでは？という点については、ノイズ自体は満遍なく広い帯域に広がっていますので、ギブス現象が顕著に現れないのではないかと思われる点
3. 周波数が高いところで減衰量が急激に大きくなるような窓関数（たとえばtukey窓…マイナーらしいですが）などが使われているかも？という点

などが思いつきました。アベレージングした結果？とかも思いましたが、FFTビン数（N）とフロアのレベル、AD9626のSNRのスペックから計算しても、そういうこともなさそうでした…」

 

ギブス現象はこの現象とは違うものだった！

しかしギブス現象は違う振る舞いでありました…。ギブス現象とは「急峻に変化する周期関数をフーリエ級数で展開し、それを有限項で打ち切ると発生する（時間信号の暴れ）[1]」とあります。ということはここで考える事象とは違うものを表しているわけですね（汗）。

ということでより詳しく「デジタル信号処理の基本」…というか奥底に立ち返って考えてみました。そのことは本稿の中盤で示していきたいと思います。

 

いぶし銀のエンジニアはさすが！

この件を知り合いの方に聞いたところ、こんなコメントもいただきました。

……窓関数を使わずに、シャープなスペクトルを得るサンプリング方法を「コヒーレント・サンプリング」と言うのですよ。たとえばFFT時間長TFFT（TFFT = N / fS）とし、

fIN ： 被測定正弦波周波数
fS ： サンプリング周波数
N = 2K ： FFTサンプル数
P ： N (= 2K)とは互いに素な数。これは被測定信号の波数（サイクル数）に比例する

としてみます。最後のPについては（互いに素な数ということもあり）奇数サイクルとなります。このようにして次の関係が成り立つようにします。

fIN / fS = P / N = P / 2K (1)

AD9626の場合、おそらく次のような設定ではないかと思いますよ。

fS = 170.00000MHz
K = 16 (N = 65536)
P = 3971
fIN = 10.3007507MHz

FFTを使ってA/Dコンバータを評価する場合、被測定信号である正弦波の位相に対して、偏りのないサンプリング（ランダムなサンプリング）をする必要があります。もしPと2Kが互いに素でないと、重複サンプリング（繰り返して同じ位相部分をサンプリングする）となり、データに無駄が生じ、かつ量子化ノイズも白色（White）にならず、一部が高調波歪に姿を変えてしまうんですよね……。

 

FFT時間長に合わせて正弦波の周波数を決めるのだ！

このコメントには「なるほど！」というところでした。FFT時間を信号の周期で「割り切れる」という視点で、信号の周波数を若干ずらして、FFT時間長（385.5059us @170Msps 64k = 65536point FFT）において、連続した（前縁と後縁で不連続が無い）信号とし、これに対してFFTを行うということなのですね！

この方は「いぶし銀」という感じの経験豊富なエンジニアなのですが、さすが！凄いなあと思ったものでした。

 

デジタル信号処理の奥底に立ち返って

FFTするということは、「被測定信号fINを、FFTとして計算処理する時間長TFFT（実際はこれまでの説明のように、A/D変換サンプリング周波数fSとFFTポイント数Nから決まる時間。TFFT = N / fS）の孤立矩形波と乗算する」と考えることができます。これを図2に示します。これは（連続時間での）フーリエ変換で周波数軸に変換したもので考えます。

窓関数を用いてFFT処理する実際の場合は、「時間TFFTの孤立矩形波と乗算する」を「時間TFFTの窓関数と乗算する」と読み替えればよいことになります。ここで被測定信号を

s(t) = A sin(2πfIN t) (2)

孤立矩形波を

b(t) = {■(1&|t|<t_FFT@0&others)┤ (3)

としてみると、この掛け算は

sig(t) = s(t)・b(t) (4)

であり、sig(t)をフーリエ変換すれば良いという考えです。

図2. ある時間長でFFTするということは被測定信号と孤立　矩形波の乗算から成り立つと考えられる

 

それぞれの波形は周波数スペクトルとしてはどうなる？

それではこのs(t), b(t), sig(t) = s(t)・b(t)は、それぞれ周波数軸上ではどのようにして（フーリエ変換として）表すことができるでしょうか。

s(t)は単一周波数なので、図3のように横軸を周波数としてみると、周波数±fINに孤立（輝線）スペクトルがあるように表すことができます。数式で表せば、

s(f) = A δ(f - fIN) + A δ(f + fIN) (5)

このうち第1項が+fINの成分（右のスペクトル）、第2項が-fINの成分（左のスペクトル）です。δはデルタ関数です。

図3. 非測定正弦波信号s(t)をフーリエ変換して周波数軸上でs(f)として表してみる

また孤立矩形波の方は、周波数軸上では図4のようになります。数式で表せば、

b(f) = 2 sin(2πf TFFT /2) / 2πf (6)

これを「sinc関数」とも呼びます。この関数は同図のように1/TFFT [Hz]の間隔でゼロをクロスします（これが以降の説明で重要）。図ではTFFT = 1secとしてあります。

 

 

図4. 孤立矩形波d(t)をフーリエ変換して周波数軸上で　　　d(f)として表してみる（TFFT = 1secとしてある）

 

時間軸上で掛け算された波形（目的の信号）は周波数スペクトルとしてはどうなる？

次にいよいよsig(t) = s(t)・b(t)ですが、これは「畳み込み」という技を使います。2つの異なる時間信号s(t), b(t)があり、これを時間軸で掛け合わせた波形sig(t) = s(t)・b(t)は、s(t), b(t)それぞれを個別にフーリエ変換して周波数軸に変換したスペクトルs(f), b(f)を畳み込みしたものになります。

ところで2つの異なる時間信号s(t), b(t)があり、s(t), b(t)それぞれを個別にフーリエ変換して周波数軸に変換したスペクトルs(f), b(f)を、周波数軸で掛け合わせたスペクトルsig(f) = s(f)・b(f)は、時間信号波形s(t), b(t)を時間軸で畳み込んだものになります。こちら（時間軸での畳み込み）はデジタル・フィルタの考え方でよく出てくるものです。

このように畳み込みとは、時間軸と周波数軸の間で、掛け算と対をなす考え方です。

このことから、s(t), b(t)を掛け合わせた周波数スペクトルというのは、図5の下側の図のように表すことができます。s(t)の周波数は5Hzとしてあります。

同図の上側の図は、右の赤い線がAδ(f - fIN)が畳み込まれたスペクトル（f = +fIN）、左の青い線が Aδ(f + fIN)が畳み込まれたスペクトル（f = -fIN）になります。

 

図5. s(t), b(t)を掛け合わせた信号sig(t)というのはそれぞれの周波数スペクトルs(f), b(f)を畳み込みしたものになる。s(t)の周波数は5Hz

畳み込みとは…

畳み込みの概念の説明をするとかなり長くなってしまいます。信号処理の参考書や参考文献[4]に記載がありますので、詳しくはこちらをご覧いただければと思います。数式で表すと

f_1 (t)⊗f_2 (t)= ∫_(-∞)^∞▒〖f_1 (x)f_2 (t-x)dx〗 (7)

となりますが、これでは何が何だか分かりませんね。イメージとして（かなり乱暴ですが）考えてみると、図6のように「ある長さのホームを超長い電車が通過するときに、『瞬間・瞬間でホーム全長の区間には何人乗っている？』を連続して観測する」と考えるようなものです。これでもわかりづらいかもしれません（汗）…。やはり[4]あたりを（特に簡単に説明してありますので）ご覧になるとよろしいかと思います。

 

図6. 畳み込み演算は、ある長さのホームを超長い電車が通過するとき、瞬間・瞬間でホーム全長の区間には何人乗っている？を連続して観測するようなもの

 

ご質問を解きほどく

さてそれでは畳み込みの概念を用いて、当初のご質問（問題）を解きほどいてみましょう。各記号の単位もあらためて表記してみます。

 

FFT時間長と被測定信号の波数（サイクル数）がぴったり合っている場合

まずは正弦波がFFT時間長で連続（正確にはサンプリングされた波形のオシリとアタマが連続）であるとき、つまり被測定信号の周波数fIN [Hz]が

fIN [Hz] = PfBIN [Hz] = P / TFFT [s] (8)

としてFFTで得られる周波数ステップである「ビン周波数」fBIN [Hz]（fBIN = 1/TFFT [s]）の整数倍（P倍）となる条件で考えてみましょう。式(1)が成り立つ条件なわけですね。確認のため式(1)を変形しておくと、

fIN [Hz] = PfS / N = PfS / 2K = P/TFFT = P fBIN (9)

となります。

この条件でのFFTのようすを連続時間として仮定してみると図5がそのまま得られるわけです。

 

次にビン周波数ステップでサンプリングされたことを考えてみる

今度これを規定のサンプリング周期（TS [s] = 1/fS）で規定の時間（TFFT [s] = N / fS）サンプリングした場合を考えてみます。サンプリングした状態でも、もともとの信号波形が周波数軸上にスペクトルとして変換されたものは連続時間の場合と同じになります。

そうすればデジタル信号処理の結果としては、図7（図5の再掲と加筆）のように、もともとの（sig(f) = s(f) ⊗ d(f)として畳み込まれた）スペクトル形状を、ビン周波数fBIN [Hz] = 1/TFFT [s]の周波数ステップで（周波数軸上で）サンプリングした●の点（離散周波数）として考えることができるわけです。これがFFTされた結果になります。

ここでさきほどの「この関数は図4のように1/TFFT [Hz]の間隔でゼロをクロスします（これが以降の説明で重要）」の説明と式(7)のとおり、fIN [Hz] = PfBINですから、sig(f)の±fIN [Hz]の周波数以外のところに存在する（連続周波数の）スペクトルは、0Hzを中心としてfBIN [Hz]間隔でサンプリングされることにより（fINはfBINステップでゼロ・クロスしているので）「全てゼロ」になります。

これにより図1で示すAD9626のデータシートのFig. 11の条件が成立していることになります。

図7.図5の信号同士が掛け合わされたスペクトルを規定の周波数間隔でサンプリングされたものとして見てみる。これが本来のFFTの結果（s(t)の周波数は5Hz）

 

FFT時間長と被測定信号の波数（サイクル数）がぴったり合っていない場合

次に被測定正弦波信号fIN [Hz]が、波形の途中で打ち切られた場合（正確にはサンプリングされた波形のオシリとアタマが連続になっていない）信号であるとき、つまり

fIN [Hz] ≠ PfBIN [Hz] = P / TFFT [s] (10)

の条件で考えてみましょう。式(10)をこんなふうに書き直してみます。

fIN [Hz] = PfBIN + d [Hz] (11)

周波数がd [Hz]だけずれているイメージです。この周波数がdだけずれたs(t)とb(t)を掛け合わせた信号sig(t)というのも、図5と同じプロセスで表すことができます。

そのようす（畳み込みの結果）が図8です。この例ではオフセット周波数d = 0.2Hzとしてあります。右の赤い線がAδ(f - PfBIN - d)が畳み込まれたスペクトル（f = +PfBIN + d）、左の青い線が Aδ(f + PfBIN + d)が畳み込まれたスペクトル（f = - PfBIN - d）になります。ここでδはデルタ関数です。

 

図8. 周波数がd = 0.2Hzだけ離れたs(t)とb(t)を掛け合わせた信号。sig(t)はさきほどから周波数dだけずれている

 

 

これをさきほどと同じようにビン周波数ステップでサンプリングしたことを考えてみる

これをさきほどと同じように、規定のサンプリング周期（TS [s] = 1/fS）で規定の時間（TFFT [s] = N / fS）サンプリングした場合として考えてみます。ここでも図9（図8の再掲と加筆）のように、もともとのスペクトル形状sig(f)をfBIN [Hz] = 1/TFFTのビン周波数ステップで（周波数軸上で）サンプリングしたものとして考えることができるわけです。

被測定信号の周波数fIN [Hz] = PfBIN + dですから、このsig(f) = s(f) ⊗d(f)として畳み込まれた信号は±dだけ外側にずれています。このようすは図9に加筆してあるとおりです。

 

図9. 図8のスペクトルを規定の周波数間隔でサンプリングしてみる（これがFFTの結果）。オフセットd = 0.2Hzにより　余計なスペクトルが出てしまう！

 

このオフセットにより、sig(f)の±fINの周波数以外のところに存在する（連続周波数の）スペクトルは、サンプリングが0Hzを中心としてfBIN間隔で行われることにより、M fBIN（M = 0～N/2）の各ポイント、つまり●の各点では「ゼロにはなりません」。

この図9はsinc関数として「電圧量」として表してみたものです。実際にスペアナなどで観測される量は「電力」です。そこで、

P [W] = V2/R（ここでR = 1Ωとして）=V2 (12)

として1Ωの抵抗に生じる電力量P [W]として計算しなおし、またこれをdBに直してプロットしたものを図10に示します。

 

図10. 図9のスペクトルを電力量Pとして、またdBにして表示してみた

もともとのご質問の答え：正弦波の途中で打ち切られた波形をFFTするとスペクトルはおかしくなる

もともとのご質問は「サンプリングが正弦波周期の途中で打ち切られる（正確にはサンプリングされた波形のオシリとアタマが連続になっていない）とおかしくならないか？」ということでしたが、答えとしては「おかしくなる」ということなのでした。

 

波数が合わないようすをMATLABでシミュレーションしてみた

正弦波のサンプリングを、① FFT時間長で連続（正確にはサンプリングされた波形のオシリとアタマが連続）であるときと、② 波形の途中で打ち切られた場合（正確にはサンプリングされた波形のオシリとアタマが連続になっていない）の差異について「シミュレーションで見てみたいなあ」と、ごにょごにょMATLABを使ってやってみました。まずはソース（mファイル）のご紹介です…（図11）。

 

% set freq parameter sf = 170E6; % sampling frequency fftpoint = 65536; ocf = 10.3E6; % output frequency (coarse setting) % change this value! phasetrunc = 0; % truncation point in degree % generate real (desired) freq number_of_wave = round(ocf * fftpoint / sf); %# of wave in FFT length, rounded period_in_point = round(fftpoint / number_of_wave); trunc_length = round(phasetrunc * period_in_point / 360) mf = number_of_wave /((fftpoint - trunc_length)/sf) %real frequency % set X axis freq = linspace(0, sf, fftpoint + 1); freq = freq(1:fftpoint); % perform FFT sampletime = [0:1/sf:(fftpoint-1)/sf]; %sample point time signal = sin(2 * pi * mf * sampletime); %sinusoidal signal spectrum = fft(signal); spectrum = 20 * log10(abs(spectrum)); %convert to dB spectrum = spectrum - max(spectrum); %normalize plot(freq, spectrum) grid on axis([0 170E6 -300 0])

図11. 波数が合わないようすをMATLABでシミュレーションしてみたmファイル

mファイル上の変数phasetrunc（黄色でハイライト）を、打ち切りたい位相量に設定すると、それを丸めた結果として周波数を設定し、FFT結果を表示します。ピークを0dB、ダイナミックレンジは300dBにしてみました。

FFTポイント数65536、fS = 170MHz、信号周波数は10.3MHz付近です（AD9626データシートFig. 11と同条件です）。mファイルを実行すると「ズレ」サンプル数とそのときの周波数が出力されます（「;」を入れてないので）。

 

シミュレーション結果の出だし

図12は正弦波が連続でFFT時間長にきちんと収まっている場合です。被測定信号の周波数はfIN = 1.030075073242188e+007 Hzになります。本来であればノイズフリーのはずですが、MATLABの計算精度によるノイズが出ていると考えられます。

図13は正弦波が最後のところが20°程度の位相で打ち切られた場合です。このとき1サンプル分の「ズレ」になります。被測定信号の周波数はfIN = 1.030090791180285e+007 Hzです。図12との差異はなんと250Hz程度です…。こんなシミュレーションはしたことがありませんが、とても興味深いものでした！

 

 

図12. 正弦波が連続でFFT時間長にきちんと収まっている場合です。被測定信号の周波数は1.030075073242188e+007 Hz		図13. 正弦波が最後のところが20°程度の位相で打ち切られた場合。被測定信号の周波数は1.030090791180285e+007 Hz

切り出し点の位相を20°ステップでずらしてシミュレーションしてみた

位相の切り出し位置をFFTでの1ポイント・ステップでずらしていき、シミュレーションしてみました。こうすると計算上、FFTでの1ポイントずれに対して20°ステップずつで、位相の不連続量が変化するように設定されます（図14～図30）。それでも割り切れないことにより、途中で10°ずれる（180°から210°へのステップが20°変化ではなく、30°変化になってしまっている）ところは出てしまっています。

下限を-120dBで切ってみました。なかなか興味深いです…。

参考文献
[1] 中村尚吾; ビギナーズデジタルフィルタ, 東京電機大学出版局
[2] 中村尚吾; ビギナーズデジタルフーリエ変換, 東京電機大学出版局
[3] 萩原将文; デジタル信号処理, 森北出版
[4] 石井聡; 合点！電子回路超入門, CQ出版社

 

図14. FFT 時間長に被測定信号がきちんと収まった場合（0°位相で切り出し。FFT ポイントとして0 ポイントのずれ）		図15. 20°位相で切り出しされた場合。FFTポイントとして1ポイントのずれ
図16. 40°位相で切り出しされた場合。FFTポイントとして2ポイントのずれ		図17. 60°位相で切り出しされた場合。FFTポイントとして3ポイントのずれ
図18. 80°位相で切り出しされた場合。FFTポイントとして4ポイントのずれ		図19. 100°位相で切り出しされた場合。FFTポイントとして5ポイントのずれ
図20. 120°位相で切り出しされた場合。FFTポイントとして6ポイントのずれ		図21. 140°位相で切り出しされた場合。FFTポイントとして7ポイントのずれ
図22. 160°位相で切り出しされた場合。FFTポイントとして8ポイントのずれ		図23. 180°位相で切り出しされた場合。FFTポイントとして9ポイントのずれ
図24. 210°位相で切り出しされた場合。FFTポイントとして10ポイントのずれ		図25. 230°位相で切り出しされた場合。FFTポイントとして11ポイントのずれ
図26. 250°位相で切り出しされた場合。FFTポイントとして12ポイントのずれ		図27 270°位相で切り出しされた場合。FFTポイントとして13ポイントのずれ
図28 290°位相で切り出しされた場合。FFTポイントとして14ポイントのずれ		図29 310°位相で切り出しされた場合。FFTポイントとして15ポイントのずれ
図30 330°位相で切り出しされた場合。FFTポイントとして16ポイントのずれ		図31 350°位相で切り出しされた場合。FFTポイントとして17ポイントのずれ

 

 

 

 

著者について

石井 聡

1963年千葉県生まれ。1985年第1級無線技術士合格。1986年東京農工大学電気工学科卒業、同年電子機器メーカ入社、長く電子回路設計業務に従事。1994年技術士（電気・電子部門）合格。2002年横浜国立大学大学院博士課程後期（電子情報工学専攻・社会人特別選抜）修了。博士（工学）。2009年アナログ・デバイセズ株式会社入社、現在に至る。2018年中小企業診断士登録。
デジタル回路（FPGAやASIC）からアナログ、高周波回路まで多...