プログラムでΔΣ変調器
by K.IIndex
ΔΣ（デルタシグマ）変調器は、アナログ信号を量子化する際に、量子化誤差をフィードバックする回路で、
そのため低域のノイズは少なくなり、量子化誤差は高域に偏る特性がある。（ノイズシェーピング）
プログラムで、ΔΣ変調器による単純な１ビットの量子化をやってみる
[top]
fftwでいろいろ実験 で、いろいろ基本的なプログラムを作ってあるので、それを使うんだけど、
今回は、FPGAとかで回路にすることを考えて、波形出力は基本的に整数値で扱う様に一部作り直す
係数を掛ける時も、出来るだけｘ２とか、ｘ４のようにシフト演算で出来るような値でやろうと思う
sin関数による正弦波生成プログラム（isingen.c）
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
//１秒分のデータを生成
#define FREQ (1000) //生成周波数
#define SAMPLE (1000000) //サンプリング周波数
#define BWIDTH 16
main( int argc, char *argv[] )
{
double d;
int f,s;
int out,t,max;
f = FREQ;
s = SAMPLE;
if (argc>1)
f = atoi(argv[1]);
if (argc>2)
s = atoi(argv[2]);
d = 2*M_PI*f/s;
max = (int)pow(2,BWIDTH-1)-1;
for (t=0; t<s ;t++) {
out = (int)(sin(t*d)*max);
printf("%d\n",out);
}
}
最初にsinで波形テーブルを作って、それを間引きして出力する。（idds.c）
生成可能な周波数fは、（サンプリング周波数S/テーブルの大きさN）の整数倍だけとなる
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
//１秒分のデータを生成
#define FREQ (1000) //生成周波数
#define SAMPLE (1000000) //サンプリング周波数
#define N (1000000) //テーブルの大きさ
#define BWIDTH 16
int sin_table[N/4+1];
int max;
void make_table(int n)
{
int x;
max = (int)pow(2,BWIDTH-1)-1;
for (x=0; x<=n/4 ;x++)
sin_table[x] = (int)(sin(2*M_PI/n*x)*max);
}
int sin_x(int x, int n)
{
int xx;
xx = x%(n/4);
if (x < n/4*1) return sin_table[xx];
else if (x < n/4*2) return sin_table[n/4-xx];
else if (x < n/4*3) return -sin_table[xx];
else return -sin_table[n/4-xx];
}
main( int argc, char *argv[] )
{
int i, d, t=0;
int f, n, s;
f = FREQ;
s = SAMPLE;
n = N;
if (argc>1)
f = atoi(argv[1]);
if (argc>2)
s = atoi(argv[2]);
n = s;
make_table(n);
d = (long)f*n/s;
for (i=0; i<s ;i++) {
printf("%d\n",sin_x(t,n));
t = (t+d)%n;
}
}
% isingen 1000 | head -1000 | gp_wav > isingen_1000.png
% idds 1000 | head -1000 | gp_wav > idds_1000.png
% isingen 1000 | fftw | gp_fft > isingen_1000_fft.png
% idds 1000 | fftw | gp_fft > idds_1000_fft.png
全然、差がなかった。。当然、整数化で切り捨てるのでdoubleで計算した時よりノイズが多くなってるけど。
sin関数直接もDDSも、同様に整数化されているので、結局全く同じ出力になっている様だ。
どちらを使っても良さそうなので、とりあえずいろいろ弄りやすいDDSの方を使おうと思う。
係数積算用のフィルタmulの整数化出力版のプログラムも作っておく（imul.c）
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
main( int argc, char *argv[] )
{
FILE *fp;
double in1,in2;
fp = fopen(argv[1],"r");
if (fp) {
while((fscanf(stdin,"%lf",&in1)!=EOF)&&(fscanf(fp,"%lf",&in2)!=EOF)) {
printf("%d\n",(int)floor(in1*in2));
}
}
else {
in2 = 1;
if (argc>1)
in2 = atof(argv[1]);
while((fscanf(stdin,"%lf",&in1)!=EOF)) {
printf("%d\n",(int)floor(in1*in2));
}
}
}
整数入力のフィルタを使うときは、これを通して整数化してから、入力すれば良い
fftw出力の表示用スクリプトは、比較しやすいようにスケール固定版を使う（gp_fftx）
#!/home/blueseed/bin/gnuplot
set terminal png
set nokey
set xlabel "f[kHz]"
set ylabel "|F(w)|"
set logscale x
set logscale y
plot [0.5:][1:1000000] '< cat -' using ($1/1000):4 with impulses linewidth 1
[top]
ΔΣ変調器は基本的に、アナログ‐デジタル変換器なので、
% idds 1000 | comp 0 | head -1000 | gp_wav > comp1000.png
% idds 1000 | comp 0 | fftw | gp_fftx > comp1000_fft.png
矩形波なので、奇数倍の高調波が含まれているが、さらに偶数倍の高調波成分も含まれている。
Dutyのズレでも偶数倍の高調波は出るけど、これは周波数によらずに同じ大きさになっている。
これは量子化誤差だと思う。全周波数領域で、均等に分布している。
#include <stdio.h>
#include <math.h>
#define BWIDTH 16
int main(int argc, char *argv[])
{
int in,out,max;
int delta=0,sigma=0;
int a1=1;
max = (int)pow(2,BWIDTH-1);
if (argc>1)
a1 = atoi(argv[1]);
while (scanf("%d", &in) != EOF) {
sigma += in - (delta ? +a1*max : -a1*max);
out = (sigma > 0) ? 1 : 0;
delta = out;
printf("%d\n",out);
}
return 0;
}
% idds 1000 | ids1 | head -1000 | gp_wav > ids1_1000.png
% idds 1000 | ids1 | fftw | gp_fftx > ids1_1000_fft_1.png
単純な２値化の場合、量子化誤差は、ほぼ全周波数に均等に分布するが、
Δシグマ変調器では、ノイズシェーピングの効果が出ていることが分かる
矩形波としての、奇数倍の高調波も少しあるようにみえる
デルタシグマ変調器の次数が大きくなれば、さらにノイズシェーピングの効果が大きくなるはず
#include <stdio.h>
#include <math.h>
#define BWIDTH 16
int main(int argc, char *argv[])
{
int in,out,max;
int delta1=0,delta2=0,sigma1=0,sigma2=0;
int a1=1,a2=1;
max = (int)pow(2,BWIDTH-1);
if (argc>1)
a1 = atoi(argv[1]);
if (argc>2)
a2 = atoi(argv[2]);
while (scanf("%d", &in) != EOF) {
sigma1 += in - (delta1 ? +a1*max : -a1*max);
sigma2 += sigma1 - (delta2 ? +a2*max : -a2*max);
out = (sigma2 > 0) ? 1 : 0;
delta1 = delta2;
delta2 = out;
printf("%d\n",out);
}
return 0;
}
% idds 1000 | ids2 | head -1000 | gp_wav > ids2_1000.png
% idds 1000 | ids2 | fftw | gp_fftx > ids2_1000_fft_1_1.png
うーん。これは発振しているのかな。。
２次のΔΣは係数をちゃんと調整しないとダメか。。。
a1=x2,x4,x8,x16で、a2=1の場合は
% idds 1000 | ids2 2 1 | fftw | gp_fftx > ids2_1000_fft_2_1.png
% idds 1000 | ids2 4 1 | fftw | gp_fftx > ids2_1000_fft_4_1.png
% idds 1000 | ids2 8 1 | fftw | gp_fftx > ids2_1000_fft_8_1.png
% idds 1000 | ids2 16 1 | fftw | gp_fftx > ids2_1000_fft_16_1.png
a1=1で、a2をx2,x4,x8,x16と大きくしてみる
% idds 1000 | ids2 1 2 | fftw | gp_fftx > ids2_1000_fft_1_2.png
% idds 1000 | ids2 1 4 | fftw | gp_fftx > ids2_1000_fft_1_4.png
% idds 1000 | ids2 1 8 | fftw | gp_fftx > ids2_1000_fft_1_8.png
% idds 1000 | ids2 1 16 | fftw | gp_fftx > ids2_1000_fft_1_16.png
a1=2で、a2をx2,x4,x8,x16と大きくしてみる
% idds 1000 | ids2 2 2 | fftw | gp_fftx > ids2_1000_fft_2_2.png
% idds 1000 | ids2 2 4 | fftw | gp_fftx > ids2_1000_fft_2_4.png
% idds 1000 | ids2 2 8 | fftw | gp_fftx > ids2_1000_fft_2_8.png
% idds 1000 | ids2 2 16 | fftw | gp_fftx > ids2_1000_fft_2_16.png
a1=4で、a2をx2,x4,x8,x16と大きくしてみる
% idds 1000 | ids2 4 2 | fftw | gp_fftx > ids2_1000_fft_4_2.png
% idds 1000 | ids2 4 4 | fftw | gp_fftx > ids2_1000_fft_4_4.png
% idds 1000 | ids2 4 8 | fftw | gp_fftx > ids2_1000_fft_4_8.png
% idds 1000 | ids2 4 16 | fftw | gp_fftx > ids2_1000_fft_4_16.png
a1＜a2にしないと、安定しないようだ。
でもa2が大き過ぎても良くない様で、a2はa1の2倍ぐらい1 a1=1の時は基本的に安定しているが、a1=2の方がちょっとノイズシェーピングの効果が大きいみたいだ
% idds 1000 | ids1 | fftw | gp_fftx > ids1_1000_fft_1.png
% idds 1000 | ids1 2 | fftw | gp_fftx > ids1_1000_fft_2.png
低域で僅かに良くなっている。とりあえず、a1=2の方がちょっとだけ良いかな2 ということで、とりあえず１次はa1=2、2次はa1=2,a2=4にしておこう
明らかに、1次より2次の方が、量子化誤差が高域に押しやられていることがわかる。
あと、2次の方がいろいろな周波数成分が生成されているようにみえる
1 といっても、あまり細かく見ていないので、もっと良い条件はあるかもしれないが。まぁ、係数の一例として。 2 高域の形が変なので、もしかすると不安定なのかもしれないが。 [top]
元の回路の積分器を遅延回路で置き換えると、こうなるので、
さらに遅延回路を２値化の前に持ってきて、積分器の遅延回路を使って省略したという感じだろうか。
これで良いのか分からないが、ループとしてみれば同じなのかな？
これに合わせて、1次のΔΣ変調器のプログラムを書き換えてみる（ids1x.c）
#include <stdio.h>
#include <math.h>
#define BWIDTH 16
int main(int argc, char *argv[])
{
int in,out,max;
int sigma=0;
int a1=1;
max = (int)pow(2,BWIDTH-1);
if (argc>1)
a1 = atoi(argv[1]);
while (scanf("%d", &in) != EOF) {
out = (sigma > 0) ? 1 : 0;
sigma += in - (out ? +a1*max : -a1*max);
printf("%d\n",out);
}
return 0;
}
sigmaの遅延回路は残ってるが、deltaの方の遅延回路は省略している。
% idds 1000 | ids1 2 | fftw | gp_fftx > ids1_1000_fft_2.png
% idds 1000 | ids1x 2 | fftw | gp_fftx > ids1x_1000_fft_2.png
１次と同様に２次のΔΣ変調器の回路は、やはり遅延回路が省略されて簡単になっている。
やはり元の回路の積分器を遅延回路で置き換えると、こうなので、
さらに遅延回路を２値化の前に持ってきて、積分器の遅延回路と共通にして省略したということか。
同様に、2次のΔΣ変調器のプログラムも書き換えてみる（ids2x.c）
#include <stdio.h>
#include <math.h>
#define BWIDTH 16
int main(int argc, char *argv[])
{
int in,out,max;
int sigma1=0,sigma2=0;
int a1=1,a2=2;
max = (int)pow(2,BWIDTH-1);
if (argc>1)
a1 = atoi(argv[1]);
if (argc>2)
a2 = atoi(argv[2]);
while (scanf("%d", &in) != EOF) {
out = (sigma2 > 0) ? 1 : 0;
sigma2 += sigma1 - (out ? +a2*max : -a2*max);
sigma1 += in - (out ? +a1*max : -a1*max);
printf("%d\n",out);
}
return 0;
}
% idds 1000 | ids2 2 4 | fftw | gp_fftx > ids2_1000_fft_2_4.png
% idds 1000 | ids2x 2 4 | fftw | gp_fftx > ids2x_1000_fft_2_4.png
やはり結果は全く同じ。ではなく、よく見ると僅かに違っている
ちょっとだけ、回路が簡単になる3
% idds 1000 | comp 0 | imul 255 > iddscmp.dat
% txt2wav -1000000 iddscmp.dat
% lame iddscmp.wav iddscmp.mp3
単純な矩形波なので、高調波を含むはっきりした音になっている
% idds 1000 | ids1x 2 | imul 255 > idds1x_2.dat
% txt2wav -1000000 idds1x_2.dat
% lame idds1x_2.wav idds1x_2.mp3
ノイズシェーピングの効果で、ほとんどノイズは感じられない
% idds 1000 | ids2x 2 4 | imul 255 > idds2x_2_4.dat
% txt2wav -1000000 idds2x_2_4.dat
% lame idds2x_2_4.wav idds2x_2_4.mp3
予想では、FFTの結果が1次ΔΣより良いので、ノイズが感じられないと思っていたが、
実際に聞いてみると、明らかに2次ΔΣの方が、ノイズを感じる。。
lameで、mp3に変換する際にリサンプリングされているのが原因かも4
Resampling: input 1000 kHz output 48 kHz
Using polyphase lowpass filter, transition band: 16452 Hz - 17032 Hz
3 この例では、省略したのは1bitの遅延回路なので、たいしたこと無いんだけど。 4 変換前のwavファイルを再生するとノイズが無いので、おそらくリサンプリングが原因と思っている。でもそれなら、1次ΔΣでノイズを感じない理由は何故だろう。 [top]
ちょっと中途半端ではあるので、感想に近いけど、
デルタシグマのノイズシェーピングは、回路が簡単なのに効果が大きくて、すごく不思議。
やはり高次のデルタシグマの方が、さらに効果が大きいけど、すぐに発散してしまうので、係数求めるのが難しい。
自分は試行錯誤で係数を求めているので、いろいろ問題はあるかもしれません。
サンプリング周波数と正弦波の周波数の関係で変わってくるかも
2次のΔΣの音にノイズが入るのは、MP3への変換時のリサンプリングが原因と思うが
高次の場合、いろいろな回路があるみたいだけど、自分には未だ理解できないので、今後の課題
max = (int)pow(2,BWIDTH-1)-1;
こうしないと、指定したビット幅を超えてしまう可能性があった。
結果は変わらないが、この文書中の計算は全てやりなおして、図も差換えた。
でも、いずれにせよオーバーフローの処理をしていないので、
MP3のリサンプリングでノイズが入るのは、オーバーフローによるピークリミッタが掛かっているのかもしれない5 5 と思ったけど、レベルを下げても同じなので、そういうわけでは無いらしい。 [top]
ΔΣ変調の部屋
ΔΣ変調について、詳しく解説されている
久々に見たらリンク切れ。凄く参考になったのになぁ。。復活させて欲しい（161227）
[top]
[プログラムの部屋に戻る]
