ExcelでDSPをしよう(DDS編)

さて、手始めにNCOを作ってみましたのですが、実際、離散領域で波形を作る場合はあのような構成にはなっていません。
・アキュムレータで扱える数値に上限がある
・sin関数の演算に時間がかかる
というのが理由です。
これらを解決した方法にDDS(Direct Digital Synthesizer)というものがあります。
プログラムで作ろうとすると面倒なので、Excelで正弦波を出力するDDSを作ってみましょう。

先に書いた欠点を克服するには
・上限を超えたらアキュムレータの数値を小さくする
・sin関数をあらかじめ計算して記憶しておく
とすればいいです。
一定数値で変更することにより周期ができます。
基本的な波形を数列/配列としてメモリに記憶し、数列/配列のインデックスをたよりに読み出せば波形が再生されます。
さらに、アキュムレータに与える数値を変えると出力の周波数が変わります。

言葉では分かりづらいと思いますが、組み立てていくと少しずつ分かっていくはずです。
ちなみにブロック図は以下です。
Zoom
NCOに比べ、簡単になった代わりに制約が増えている事が分かります。
本来ならまだ制約があるのですが、そこまで難しくする必要は無いので無視します。

  まずは解決法通り、sin関数をあらかじめ計算し記憶する、という工程を行います。
分かりにくくなるのであらかじめサンプル数nの最大値nmaxを決めておきます。
一般的には2^n-1の数値にします。(0から始まる2進数に基づくのでキリが良い)
今回はソフトウエア処理なので2^n-1の個数にしなくてもいいのですが、とりあえず定石に従い
nmax=1023
としましょう。
Zoom ・n:サンプル番号(配列のインデックス)
・sin[n]:sin(n)の結果
を表に加えます。
ただし、
・nは0〜1023(全部で1024個)
・sin[n]=sin(2πn/(nmax+1))=sin(2πn/1024)
です。
単位円1周である2π radを1024等分したという事です。
nとsin[n]をそれぞれ1024個分生成すれば、sin関数の配列が完成します。
このように入力と出力の関係を列挙した物を「一覧表、ルックアップテーブル(LUT)」と言います。
次にnのアキュムレータを作ります。
先ほどのLUTから離れた場所にパラメータの
・fs:サンプリング周波数 [Hz]
・a:アキュムレータへの加算値
を用意しておきます。

続いてNCOと同様に
・n:サンプル番号
・t:時間 [s]
・Vo:出力 [V]
を入力しておきます。
今回は振幅のパラメータAを省略します。
Zoom nの初期値は0です。
さらにnの値を条件に従って変更するのでIF()関数を使います。
変更条件は
n>nmax
と書けます。等号が無い事に注意しましょう。
ですが、Excelの性質上自分自身のセルを参照できないので
条件:(一つ前のn)+a>nmax
真ならば:(一つ前のn)+a-nmax
偽ならば:(一つ前のn)+a
とすれば良い事になります。
これらをExcelの書式で書きましょう。
tのアキュムレータも作ります。
これはNCOと同様に漸化式は
t[0]=0 t[n]=1/fs+t[n+1]
です。
アキュムレータではなく、
t[n]=n/fs
としてしまうとtはリセットされてしまうので使えません。
今回、tはグラフの横軸としてだけ使います。
Zoom Voである「LUTからの読み出し」を行います。
これはExcelにあるVLOOKUP()関数を使います。
この関数は指定範囲内の1列目をインデックスとし、その右側にあるセルの内容を返します。
よって、関数に与える指定パラメータは
・検索値:nのアキュムレータの値
・範囲:LUTの全範囲
・列番号:2 (1はインデックス、2はsin[n]の値)
・検索方法:FALSE (完全一致)
となります。
範囲は固定なので、絶対参照にしておきます。
a=1 fs=1000
を入力しましょう。
そして、n、t、Voの連続データを作ります。
nがリセットされる事を確認したいので、総数が1025個以上になるように作ります。
nはaずつ加算されてゆき、1023の次が0になりましたか?
Zoom 結果を確認しましょう。
tとVoをすべて選択し前回同様、散布図(直線)を作ります。
もちろん正弦波が出ます。
この時、グラフ(もしくは表)を見ると
1/f=nmax/fs
になっている事が分かります。
このままではLUTの値がそのまま出て、横軸がtに変わっただけです。
Zoom では、aを変えましょう。
a=2
に変更します。
グラフには2周期分の波形が出ましたか?
今度は自分の好きな数値を入れてください。
その数値分の周期が表示されます。
すなわち
・f=fs・a/(nmax+1)
・fmin=fs/(nmax+1) (fmin:最小周波数、周波数分解能)
となります。
当然ながら、離散信号なのでaを上げまくると波形は崩れます。

上式より、fからaを設定する事もできます。
ですが、インデックスは必ず正の整数なので、その場合はaを整数に直さないといけません。
整数に直す事により、実際の出力周波数に丸め誤差が発生することになります。

周波数分解能を上げるためにはLUTの長さ(nmax)を大きくすれば良いというのは自明です。
ですが、どのように実装しても(もちろんExcelでも)LUTはメモリに格納されます。
分解能を良くする事は、メモリを多く消費する事と同じなのです。

NCOの時、位相差について省略しました。
NCOに位相差φを付けるなら
Vo=A・sin(2πft+φ)
となるのは御存じの通りです。
ではDDSではどうでしょう。
NCOと同じように、周波数の項、すなわちnの初期値に規格化位相差Φを足せばいいだけです。
ただ、角度の単位はラジアンになってしまうので分かりにくいです。
位相差φ、一周2πに対する割合pを使い
φ=2π/p
と書けば分かりやすいはずです。
さらにφは2π/nmaxで規格化されていないといけないので
Φ=nmax/2π・φ=nmax・p
さらに先ほどの通りインデックスは整数なので
Φ=int(nmax・p)
となります。
よってブロック図は
Zoom
となります。

fs、aのパラメータ付近に
p:位相差の割合
を加えます。
整数変換の関数はExcelでも
INT()
です。
アキュムレータを改造するのではなく、初期値の所へ
Φ=int(1024・p)
を入れないといけません。
=IF(〜)+Φ
としてしまうと周波数が変わることになります。
Zoom 試しに
p=0.25
と入れてみましょう。
位相差π/2が与えられ、波形は余弦波になるはずです。
ただし、pの最大値は1未満です。
それ以上を入れると初期値がエラーになります。
しかし、IF()関数のおかげで、初めの方の値がエラーになるだけでアキュムレータは動作しつづけます。
余談:
初期値のエラー対策を行わなくても、Excelではエラーを回避できました。
ですが、実際のプログラムではポインタ/アドレスのオーバーランとなり、動作が停止します。
パソコンでのプログラミングでよく目にする「インデックスが配列の境界外です」が起きます。

おそらく、正弦波に限定せずとも、周期を持つ任意波形までも生成できる事に気が付いているはずです。
これもDDSの強みの一つです。
n=0〜nmaxを1周期、1を最大値とした波形数列をLUTに書き込めば良いですね。
ここでは書かないので、各自やってみましょう。

戻る
inserted by FC2 system