fftwでいろいろ実験

Index

• 概要
• fftwインストール
  • とにかくインストールしてみる
  • ちゃんとインストールする
  • fftwのプログラム
• gnuplotインストール
  • gnuplotのスクリプト
  • スクリプトで標準入力
  • 表示用スクリプトいろいろ
  • 位相アンラッピング
• 信号発生
  • 正弦波発生
  • DDS
  • 矩形波
  • Dutyを変えてみる
  • 三角波
  • ノコギリ波
  • ノイズ
  • PN符号
  • PN符号（１ビット）
• 信号の組み合わせ
  • 足し算
  • 掛け算
  • 比較
  • 最小値・最大値
• 音声データへの変換ツール
  • txt2wav
  • ステレオにする（170407追記）
  • dumpwave
  • lame
  • sox
• 音で確認
  • DDSの音
  • 三角波の音
  • 三角波の音（改）
  • 矩形波の音
  • ノコギリ波の音
  • ノイズの音
  • PN符号の音
  • DDSの音（サンプリング1MHz）
• ２値化
  • 単純な２値化
  • PWM
  • PWMの周波数を上げる（161223追記）
  • RPWM
  • １次ΔΣ
• 擬似ホワイトノイズを作る
  • 正弦波の加算プログラム
  • サンプリング周波数1kHz（失敗）
  • 位相特性用には使える
  • サンプリング周波数1kHz（改）
  • サンプリング周波数48kHz
  • 指定範囲だけ生成
• まとめ
  • 参考

概要

• さくらインターネットのスタンダードプランで、Webサーバを運用中
• 現在、自宅にはUnixマシンが無い。

• 高速フーリエ変換用ライブラリのfftwが使いたくなったので、さくらのサーバにインストールしてみた
• ここはFreeBSDだけど、まぁ、普通のLinuxでも同じだと思う

• fftwは、サンプル数を自由に選べて、とにかく簡単に周波数成分を調べられる
• ということで、fftwを使って、いろいろ遊んでみました。
• 主に、デジタルフィルタの勉強用かな。

fftwインストール

とにかくインストールしてみる

• fftwのダウンロードページから、ダウンロードして
• 普通の環境なら、以下でインストールできる

  % gunzip -c fftw-3.3.4.tar.gz | tar -xvf -
  % cd fftw-3.3.4
  % ./configure
  % make
  % make install
  

• 自分の環境では、ルート権限が無いので、インストールに失敗する
• でも、ライブラリ自体は出来てるので、

• 以下のファイルを、カレントディレクトリにコピーする
• .libs/libfftw3.a
• api/fftw3.h

• 一応、コンパイル出来たんだけど。。

$ gcc fftw.c -o fftw libfftw3.a -lm

• ちょっとカッコ悪い。。。

ちゃんとインストールする

• configureの時、prefixでユーザのHOMEを指定することで、
• その直下に、インストールされるみたいだ

  % ./configure --prefix=/home/username
  % make
  % make install
  

• /includeと、/libにインストールされた

• コンパイル出来ました！

% gcc fftw.c -o fftw -lm -lfftw3 -L/home/username/lib -I/home/username/include

• この方が、気分が良い。やってることは同じなんだけどね。。

• コンパイルしたfftwは、pathの通っているディレクトリ、例えばHomeの/binとかに入れておく。

• 管理者権限があれば、--prefixとかいらないので、念のため。。

fftwのプログラム

• fftwを使って、fftを行う簡単なフィルタプログラムを組んでみた
• 入力データとして、改行コード区切りの数値データを用意すれば良い
• サンプル数は、任意で構わない（行数がサンプル数となる）

  #include<stdio.h>
  #include<stdlib.h>
  #include<math.h>
  #include "fftw3.h"
  
  #define SIZE 1000000000
  
  main( int argc, char *argv[])
  
  {
          FILE *fp;
          int i;
          int N;
          double *Y; 
          double re,im,mag,ang;
          fftw_complex *in, *out;
          fftw_plan p;
  
          // データ読込み用のバッファを確保
          if (!(Y=(double *)malloc(SIZE))) exit(-1);
  
          // データ読込み（データ個数を調べる）
          for (N=0; scanf("%lf",&Y[N])!=EOF && N<SIZE ;N++) ;
          
          //N個のFFTW用の複素数配列を確保する
          in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);
          out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*N);
          
          // 複素数配列に入力データをセット
          for(i=0; i<N; i++) {            // t=0〜1secのデータ
                  in[i][0]=Y[i];          // 実数部
                  in[i][1]=0.0;           // 虚数部
          }
  
          // １次元のフーリエ変換を実行
          p = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
          fftw_execute(p); 
          
  
          // 結果を出力する
          for(i=0; i<N; i++) {            // f=0〜N（サンプリング周波数）のデータ
                  re = out[i][0];         // 複素数の実数部
                  im = out[i][1];         // 複素数の虚数部
                  mag = sqrt(re*re + im*im);      // 大きさを計算
                  ang = atan2(im,re);             // 位相角を計算
                  printf("%d \t %.5lf \t %.5lf \t %.5lf \t %.5lf \n", i,re,im,mag,ang);
          }
  
          // 後始末（使用した配列等を廃棄）
          fftw_destroy_plan(p);
          fftw_free(in);
          fftw_free(out);
  
          free(Y);
  }
  

• 位相角の計算式がimとreが逆になっていたので修正しました。（160920）

• このプログラムは、N個のデータ（改行コードで区切る）を読み込んで、
• 周波数０〜Nの範囲のFFT結果を出力する
  • 周波数、実数部、虚数部、絶対値、位相角の順で出力される

• fftというと、通常サンプル数が２の倍数である必要があったりするが、
• fftwは、２の倍数じゃなくても大丈夫なアルゴリズムらしい

gnuplotインストール

• fftwの結果を見るために、gnuplotもインストールする

• こちらを参考に
• gnuplotから、ダウンロード

• 自分のHOMEディレクトリにインストール

% gunzip -c gnuplot-4.6.5.tar.gz | tar -xvf -
% cd gnuplot-4.6.5
% ./configure --prefix=/home/username
% make
% make install

gnuplotのスクリプト

• gnuplot用にスクリプトファイルを作成する

set terminal png
set nokey
set xlabel "f[kHz]"
set ylabel "|F(w)|"
set logscale x
set logscale y
plot [:][:] 'fftw.dat' using ($1/1000):4 with impulses linewidth 1

• fftw.datファイルを入力として、
• png形式で出力
• X,Y軸はLogスケールで
• X軸は、１項め/1000（$を付ければ値を計算して出力できる）
• Y軸は、４項めをライン幅１のインパルスで出力

• つまりX軸に周波数（kHz）、Y軸にFFT結果の絶対値を出力する

• スクリプトの実行は、スクリプト名を gp_outとすると

$ gnuplot gp_out > hoge.png

スクリプトで標準入力

• gnuplotはファイル名を'-'にすると標準入力に出来るが、スクリプトファイルの場合はうまくいかない。
• でも、 こちらで、うまくやっていたので、真似させてもらう¹。

• fftw出力の表示用のスクリプト（gp_fft）

#!/home/username/bin/gnuplot

set terminal png
set nokey
set xlabel "f[kHz]"
set ylabel "|F(w)|"
set logscale x
set logscale y
plot [:][:] '< cat -' using ($1/1000):4 with impulses linewidth 1

• これで、中間ファイル無しで、直接プロットファイルを作成できる

$ singen 1000 | fftw | gp_fft > hoge.png

表示用スクリプトいろいろ

• gnuplot表示用スクリプトは、いくつか作っておいて使い分ける

• fftw出力の表示用のスクリプト（周波数リニアスケール版）（gp_fftl）

#!/home/username/bin/gnuplot

set terminal png
set nokey
set xlabel "f[kHz]"
set ylabel "|F(w)|"
set logscale y
plot [:][:] '< cat -' using ($1/1000):4 with impulses linewidth 1

• fftw出力の表示用のスクリプト（ゲイン・周波数リニアスケール版）（gp_fftll）

#!/home/username/bin/gnuplot

set terminal png
set nokey
set xlabel "f[kHz]"
set ylabel "|F(w)|"
plot [:][:] '< cat -' using ($1/1000):4 with impulses linewidth 1

• fftw位相表示用のスクリプト（gp_ang）

#!/home/username/bin/gnuplot

set terminal png
set nokey
set xlabel "f[kHz]"
set ylabel "Angle"
set logscale x
plot [:][:] '< cat -' using ($1/1000):5 with line linewidth 2

• 波形表示用のスクリプト（gp_wav）

#!/home/username/bin/gnuplot

set terminal png
set nokey
plot [:][:] '< cat -' with line linewidth 5

• 波形表示用のスクリプト（pointのみ表示）（gp_wavp）

#!/home/username/bin/gnuplot

set terminal png
set nokey
set style line 1 pointtype 3
plot [:][:] '< cat -' linestyle 1

• 波形表示用のスクリプト（impulse表示）（gp_pulse）

#!/home/username/bin/gnuplot

set terminal png
set nokey
plot [:][:] '< cat -' with impulses linewidth 5

位相アンラッピング

• 位相を表示すると、-π〜πの範囲に入らない波形の場合、不連続なグラフになってしまう。
• これだと、見にくいので、繋げて表示するルーチンを作成（fftw_unwrap.c）

  #include <stdio.h>
  #include <math.h>
  
  double last;
  
  double unwrap(double in)
  {
          double div;
  
          while (in-last < -2*M_PI)
              in += 2*M_PI;
          while (in-last > 2*M_PI)
              in -= 2*M_PI;
  
          div = in - last;
          if (div < -M_PI)
              in += 2*M_PI;
          else if (div > M_PI)
              in -= 2*M_PI;
          last = in;
  
          return in;
  }
  
  int     main(int argc, char *argv[])
  {
          int i;
          double re,im,mag,ang;
  
          while (scanf("%d %lf %lf %lf %lf",&i,&re,&im,&mag,&ang) != EOF) {
            printf("%d \t %.5lf \t %.5lf \t %.5lf \t %.5lf\n",
                  i,re,im,mag,unwrap(ang));
          }
          return 0;
  }
  

• とりあえず、fftwの出力の位相成分をunwrapするようにした
• 使い方は、こんな感じで

  % white0 1000 | lpf 0.5 | fftw | fftw_unwrap | gp_ang > lpf_ang.png
  

• 改良の余地はありそうだけど、一応使えるかな

• 雑音を含んでいたり、サンプリングが荒い場合は、位相ジャンプが起きてしまうので注意
• こちらを参考にさせてもらいました

信号発生

• 基本的な信号をプログラムで発生する
• 早速、fftwのプログラムを試してみよう

正弦波発生

• 信号源として、正弦波を発生させるプログラム（singen.c）
• サンプリング周波数1MHzで、１秒分のデータ（1000000point）

  #include<stdio.h>
  #include<math.h>
  
  #define SAMPLE 1000000
  #define MAX      32767
  #define POINT  1000000
  
  main( int argc, char *argv[] )
  {
          double d,f;
          double out;
          int t;
  
          f = atoi(argv[1]);
          d = 2*M_PI*f/SAMPLE;
  
          for (t=0; t<POINT ;t++) {
            out = sin(t*d);
            printf("%d\n",(int)(out*MAX/2));
          }
  }
  

• singenの出力を最初の5000pointだけgnuplotで出力してみる

% singen 1000 | head -5000 | gp_wav > singen.png

• 1kHzのsingen出力を、fftwを通してgnuplotで出力

% singen 1000 | fftw | gp_fft > singen_fft.png

• 1kHz以外に、奇数倍の高調波があることがわかる

• 試しに、データを一個だけ削ってみる

% singen 1000 | head -999999 | fftw | gp_fft > singenx_fft.png

• 矩形窓の影響らしきものが出てくる

• 実際のサンプリングでは、周波数をピッタリ合わせることは出来ないので、
• 窓の影響が、かなり大きいことが良く分かる。

DDS

• 別の方法として、テーブルを使った正弦波の生成を考えてみる

• まず、１周期分の正弦波のデータをテーブルにしておく
• データを間引きして出力することで周波数を変えることが出来る

• テーブルの大きさをN、サンプリング周波数をSとして、
• テーブルを読み出す間隔dを変更することで、生成される周波数fが変化する

  f = S/N * d
  
  d = N/S * f
  
  N = d/f * S
  

• 生成可能な周波数fは、（サンプリング周波数S/テーブルの大きさN）の整数倍となる
• 例えば、 サンプリング周波数S=1MHz、N=4000の場合、250Hzの整数倍の周波数の生成が可能になる

• 1/4周期分のテーブルを用意して、象限毎にデータを出力すれば良い

• 実際には、0〜N/4で、N/4+1個のテーブルを用意して、以下のようにプログラム（dds.c）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

//１秒分のデータを生成
#define FREQ   (1000)           //生成周波数
#define SAMPLE (1000000)        //サンプリング周波数
#define N      (1000*4)         //サンプル数

double sin_table[N/4+1];

void make_table(void)
{
        int x;

        for (x=0; x<=N/4 ;x++)
          sin_table[x] = sin(2*M_PI/N*x);
}

double sin_x(int x)
{
        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;
        double f, s;

        f = FREQ;
        s = SAMPLE;
        make_table();
        if (argc>1)
          f = atof(argv[1]);
        if (argc>2)
          s = atof(argv[2]);

        d = f*N/s;

        for (i=0; i<s ;i++) {
                printf("%f\n",sin_x(t));
                t = (t+d)%N;
        }
}

• 一応、最初の5000ポイントを出力して、波形を確認する

% dds 1000 | head -5000 | gp_wav > dds.png

• fftを掛けてみる

% dds 1000 | fftw | gp_fft > dds_fft.png

• 奇数倍の高調波はあるが、単純にsin関数を使った場合より高調波が少なくなっている
• おそらく、ddsの方は完全に周期的なデータなので、計算上誤差が小さくなるんだろうと思う

• これなら、実験には十分使えそうだ。
• とりあえず正弦波の生成には、ddsのルーチンを使うことにする

矩形波

• 矩形波の生成は、DDSのプログラムをちょっと弄っただけで簡単に作れる（square.c）
• 象限毎に、1, 1, -1, -1 と変化させるだけ

  #include<stdio.h>
  #include<stdlib.h>
  #include<math.h>
  
  //１秒分のデータを生成
  #define FREQ   (1000)           //生成周波数
  #define SAMPLE (1000000)        //サンプリング周波数
  #define N      (1000*4)         //サンプル数
  
  int square_x(int x)
  {
          int xx;
  
          if      (x < N/4*1) return 1;
          else if (x < N/4*2) return 1;
          else if (x < N/4*3) return -1;
          else return -1;
  }
  
  main( int argc, char *argv[] )
  {
          int i, d, t=0;
          double f, s;
  
          f = FREQ;
          s = SAMPLE;
          if (argc>1)
            f = atof(argv[1]);
          if (argc>2)
            s = atof(argv[2]);
  
          d = f*N/s;
  
          for (i=0; i<s ;i++) {
                  printf("%d\n",square_x(t));
                  t = (t+d)%N;
          }
  }
  

• 矩形波の出力を最初の5000pointだけgnuplotで出力してみる

% square 1000 | head -5000 | gp_wav > square.png

• 矩形波のFFTを掛けてみる

% square 1000 | fftw | gp_fft > square_fft.png

• 奇数倍の高調波のみ含まれている

Dutyを変えてみる

• Dutyを変えたらどうなるか、とりあえず矩形波の第３象限も１にして、Dutyを変えてみる（150608追記）
• 最初の5000pointだけgnuplotで出力してみる

  % squarex 1000 | head -5000 | gp_wav > squarex.png
  

• Dutyを変えた矩形波のFFT

% squarex 1000 | fftw | gp_fft > squarex_fft.png

• 矩形波のDutyが50%の時は、奇数倍の高調波のみ含まれているが、
  • Dutyが変わると、偶数倍の高調波も含まれる様になることが分かる

三角波

• 三角波の生成も、DDSのテーブルを変えるだけなので簡単！（triangle.c）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

//１秒分のデータを生成
#define FREQ   (1000)           //生成周波数
#define SAMPLE (1000000)        //サンプリング周波数
#define N      (1000*4)         //サンプル数

double tri_table[N/4+1];

void make_table(void)
{
        int x;

        for (x=0; x<=N/4 ;x++)
          tri_table[x] = 4.0/N*x;
}

double tri_x(int x)
{
        int xx;

        xx = x%(N/4);
        if      (x < N/4*1) return tri_table[xx];
        else if (x < N/4*2) return tri_table[N/4-xx];
        else if (x < N/4*3) return -tri_table[xx];
        else return -tri_table[N/4-xx];
}

main( int argc, char *argv[] )
{
        int i, d, t=0;
        double f, s;

        f = FREQ;
        s = SAMPLE;
        make_table();
        if (argc>1)
          f = atof(argv[1]);
        if (argc>2)
          s = atof(argv[2]);

        d = f*N/s;

        for (i=0; i<s ;i++) {
                printf("%f\n",tri_x(t));
                t = (t+d)%N;
        }
}

• 三角波の出力を最初の5000pointだけgnuplotで出力してみる

% triangle 1000 | head -5000 | gp_wav > triangle.png

• 三角波のFFTを掛けてみる

% triangle 1000 | fftw | gp_fft > triangle_fft.png

• パッと見た感じでは、矩形波と、ほとんど区別がつかない
• でも奇数倍の高調波のみ含まれているが、周波数の高い成分が急激に少なく²なっている

ノコギリ波

• ノコギリ波の生成も、テーブルを２つに分けるようにして生成できる（saw.c）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

//１秒分のデータを生成
#define FREQ   (1000)           //生成周波数
#define SAMPLE (1000000)        //サンプリング周波数
#define N      (1000*4)         //サンプル数

double saw_table[N/2+1];

void make_table(void)
{
        int x;

        for (x=0; x<=N/2 ;x++)
          saw_table[x] = 2.0/N*x;
}

double saw_x(int x)
{
        int xx;

        xx = x%(N/2);
        if      (x < N/2*1) return saw_table[xx];
        else return -saw_table[N/2-xx];
}

main( int argc, char *argv[] )
{
        int i, d, t=0;
        double f, s;

        f = FREQ;
        s = SAMPLE;
        make_table();
        if (argc>1)
          f = atof(argv[1]);
        if (argc>2)
          s = atof(argv[2]);

        d = f*N/s;

        for (i=0; i<s ;i++) {
                printf("%f\n",saw_x(t));
                t = (t+d)%N;
        }
}

• ノコギリ波の出力を最初の5000pointだけgnuplotで出力してみる

% saw 1000 | head -5000 | gp_wav > saw.png

• ノコギリ波のFFTを掛けてみる

% saw 1000 | fftw | gp_fft > saw_fft.png

• 三角波と違って、奇数倍だけでなく、偶数倍の高調波も含まれていることが分かる

ノイズ

• rand関数で、ノイズ出力するプログラムを書いてみた（rand.c）

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

//１秒分のデータを生成
#define SAMPLE (1000000)        //サンプリング周波数

double Uniform( void ){
    return ((double)rand()+1.0)/((double)RAND_MAX+2.0);
}

main( int argc, char *argv[] )
{
        int i;

        srand(777);
        for (i=0; i<SAMPLE ;i++) {
                printf("%f\n",Uniform());
        }
}

• rand関数によるノイズ出力を最初の100pointだけgnuplotで出力してみる

% rand | head -100 | gp_wav > rand.png

• randによるノイズにFFTを掛けてみる

% rand | fftw | gp_fft > rand_fft.png

• 右側は周波数軸を、リニアにしてみた

• randによる乱数は偏りがあるらしいが、FFTでは全然分からないなぁ。
• とりあえず、簡易的な実験には十分問題なさそうだ

PN符号

• PN符号は、シフトレジスタと加算だけで作成できる擬似乱数
• ハードウェアで、乱数を生成する時によく使われる

• 16bitのPN符号によるノイズ生成プログラム（pngen.c）
• 周期が最大になるように、タップを設定する必要がある →擬似ノイズジェネレータ

  #include<stdio.h>
  #include<math.h>
  
  //１秒分のデータを生成
  #define SAMPLE (1000000)        //サンプリング周波数
  
  int stage = 1023;
  
  int pn16(data,tap)
    int data,tap;
  {
          int mask,line,i;
  
          data = data & stage;
          mask = data & tap;
  
          line = 0;
          for (i=0; i<24 ;i++) {
            line += mask & 0x01;
            mask = mask>>1;
          }
  
          data = data<<1;
          if (line & 0x01) {
            data = data|1;
          }
  
          return(data&stage);
  }
  
  main( int argc, char *argv[] )
  {
          int     pn=0x3ff;       // 初期値（全ビット１）
          int     i;
  
          for (i=0; i<SAMPLE ;i++) {
            pn = pn16(pn,0x204);
            printf("%d\n",pn&0x3ff);      // 乱数出力
          }
  
  }
  

• PN符号によるノイズ（最初の100個のデータを出力）

% pngen | head -100 | gp_wav > pngen.png

• 基本、シフトレジスタだから、こんな波形になるのか。。

• pngenによるノイズにFFTを掛けてみる

% pngen | fftw | gp_fft > pngen_fft.png

• PN符号によるノイズはFFTを掛けると、低周波では、けっこうクセがありそうだが、
• 周波数をリニアスケールにすると、まぁそれなりで、用途によっては問題ないかも

PN符号（１ビット）

• PN符号はシフトレジスタをフィードバックして作っているので、１ビットだけ取り出す方が良いのかも。
• ということで、bit0だけ出力してみよう

• bit0だけ取り出すフィルタを作る（bit0.c）

#include<stdio.h>
#include<math.h>

main( int argc, char *argv[] )
{
        unsigned int in;

        while (scanf("%d", &in) != EOF) {
          printf("%d\n",in&0x01);
        }
}

• pngen出力を、このフィルタbit0を通して、最初の100個のデータ出力

% pngen | bit0 | head -100 | gp_wav > pngen0.png

• FFTを掛けてみると、意外と特性に大きな変化は無い
• でも、全体的には、ほぼフラットな特性になっている

  % pngen | bit0 | fftw | gp_fft > pngen0_fft.png
  

• 矩形波だと、高調波が大きいので、ノイズとしては却って良いのか？
• 1bitだけで使う方が、良いのかもしれないなぁ。

信号の組み合わせ

• ２つの信号の組み合わせ、例えば足し算とか、掛け算を考える

• 本当は、別々に作成した正弦波を足し算するようなフィルタを作れたら良いんだけど
• シェルで、２つの標準入力を動的に指定することは出来ない。

• まぁ、これが出来たら、多分マルチスレッドになっちゃう気がするので、無理なんだろう
• 片方をスタティックなファイルにして、標準入力と演算するか、
• 或いはプログラムで生成したデータを、標準入力と組み合わせるしかないかな〜

足し算

• 中間ファイルが必要になるけど、ファイルの内容と標準入力を加算するフィルタを作ってみる（add.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("%f\n",in1+in2);
            }
          }
          else {
            in2 = atof(argv[1]);
            while((fscanf(stdin,"%lf",&in1)!=EOF)) {
              printf("%f\n",in1+in2);
            }
          }
  }
  

• まず、1kHzの正弦波データを作っておいて、
• それに、2kHzの正弦波データを加算した波形の最初の5000ポイントだけ出力する

  % dds 1000 > dds1000.dat
  % dds 2000 | add dds1000.dat | head -5000 | gp_wav > add_1000_2000.png
  

• fftを掛けてみる

% dds 2000 | add dds1000.dat | fftw | gp_fft > add_1000_2000_fft.png

• 和の場合は、単純に２つの周波数が合成されているだけだ

掛け算

• 加算と全く同じだけど、ファイルの内容と標準入力を積算するフィルタを作ってみる（mul.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("%f\n",in1*in2);
            }
          }
          else {
            in2 = atof(argv[1]);
            while((fscanf(stdin,"%lf",&in1)!=EOF)) {
              printf("%f\n",in1*in2);
            }
          }
  }
  

• で、1kHzの正弦波データを作っておいて、
• それに、2kHzの正弦波データを積算した波形の最初の5000ポイントだけ出力して
• さらにfftを掛けてみる

  % dds 1000 > dds1000.dat
  % dds 2000 | mul dds1000.dat | head -5000 | gp_wav > mul_1000_2000.png
  % dds 2000 | mul dds1000.dat | fftw | gp_fft > mul_1000_2000_fft.png
  

• 1kHzと2kHzの、和である3kHzと、差分の1kHzが出力されている

• なんか分かり難いので、1kHzと1.25kHzで、同様に積算してみた

% dds 1250 | mul dds1000.dat | head -5000 | gp_wav > mul_1000_1250.png
% dds 1250 | mul dds1000.dat | fftw | gp_fft > mul_1000_1250_fft.png

• 1kHzと1.25kHzの、和である2.25kHzと、差分の0.25kHzが出力される

• 積算では、２つの周波数の和と差分が出てくるので、
• 掛け算によって、信号が周波数変換されていることが分かる

比較

• 同様に、ファイルの内容と標準入力を比較するフィルタを作ってみる（comp.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",in1>in2);
            }
          }
          else {
            in2 = atof(argv[1]);
            while((fscanf(stdin,"%lf",&in1)!=EOF)) {
              printf("%d\n",in1>in2);
            }
          }
  }
  

• 試しに、三角波に対して、0.5の閾値で比較してみた
• 三角波を最初の5000pointだけ出力して比較する

  % triangle 1000 | head -5000 | comp 0.5 | gp_wav > comp05.png
  

最小値・最大値

• ２つの信号の組み合わせじゃないんだけど、
• 最小値、最大値を求めるフィルタのプログラム（minmax.c）

  #include<stdio.h>
  #include<stdlib.h>
  #include<float.h>
  #include<math.h>
  
  main( int argc, char *argv[] )
  {
          FILE *fp;
          double in, min=DBL_MAX, max=-DBL_MAX;
  
          while((fscanf(stdin,"%lf",&in)!=EOF)) {
            if (min>in) min = in;
            if (max<in) max = in;
          }
          printf("%f, %f\n",min,max);
  }
  

• 使い方は、標準入力から信号を入れるだけ

% triangle 1000 | minmax
-0.997000, 0.997000

• 最小値、最大値が表示される

• データ範囲の確認用として作成した

• だんだんと道具が揃ってきたかな

音声データへの変換ツール

• 波形を見ていても、なかなか違いが分からない。
• 音で聞けば、分かりやすいかもということで、サウンドファイルに変換する方法を考える

txt2wav

• コマンドラインからwaveファイルを作る方法の、 txt2wav.cというプログラムを使わせてもらった。
• でも、コンパイルしたプログラムで変換すると、ヘッダのエラーになる

• 調べてみると、生成されたサウンドファイルのヘッダが44byteより大きくなってエラーになっている
• どうも、64bit環境でコンパイルしたので、ヘッダが大きくなりすぎてしまう様だ

• 試しに、32bitでコンパイルしてみたら変換できるようになった

# gcc -m32 txt2wav.c -o txt2wav -lm

• 試しに、1kHzの正弦波のwavファイルを生成。サンプリング周波数は50kHz
• dds出力は-1〜1なので、+1してx128して切り捨てることで、0〜255になるように調整

  % dds 1000 50000 | add 1 | mul 128 | int > dds1000x.dat
  % txt2wav -50000 dds1000x.dat
  

• これで、dds1000x.wavが出来る →聞いてみる
• 再生時に別ページが開くのが鬱陶しいので、 bgsoundを使う方法もあるんだけど、これはIEしか使えない

• どんなブラウザでも再生できるようにするには、 audio要素を使うべきだろうか
• でも、IEだとWAVは再生できないので、MP3にしなきゃダメみたい。。。あぁ面倒だ。。。

ステレオにする（170407追記）

• 作成したサウンドファイルはモノラルなので、ステレオにしてみる。
• といっても、単に同じデータを並べるだけなので、本当のステレオではないんだけど（stereo.c）

  #include<stdio.h>
  #include<stdlib.h>
  #include<math.h>
  
  main( int argc, char *argv[] )
  {
          double in;
  
          while((fscanf(stdin,"%lf",&in)!=EOF)) {
              printf("%d\t%d\n",(int)floor(in),(int)floor(in));
          }
  }
  

• 正弦波データを並べたデータを作って、txt2wavに掛ける
• txt2wavは、ステレオデータに対応しているので簡単

  % dds 1000 44100 | mul 32767 | stereo > dds1000s.dat
  % txt2wav -44100 dds1000s.dat
  

dumpwave

• 奥村晴彦さんの WAVE_file_formatのページから、dumpwave.cというプログラムを使わせてもらう

% dumpwave dds1000x.wav > dds1000x.txt

• これで、逆にwavファイルから、データファイルを作って、
• データ処理の実験をすることも可能になった

lame

• Linuxで動くフリーのMP3エンコーダに、 lame（Lame Ain't an MP3 Encoder）というのがあるらしい
• 早速ダウンロードして、インストールしてみた（自分の環境はFreeBSDだけど）

  % gunzip -c lame-3.99.5.tar.gz | tar -xvf -
  % cd lame-3.99.5
  % ./configure --prefix=/home/username
  % make
  % make install
  

• 簡単だった。。

• wavファイルを、MP3に変換してみる

% lame dds1000x.wav dds1000x.mp3

• Audio要素は、ちょっと大げさだけど、まぁいいか。。

• コントロールの表示方法は変えられるかもしれないけど、まぁそれが目的じゃないから。

sox

• sox（Sound eXchange）というのもあるらしい。（140902追記）
• 結構、いろんなフォーマットに対応している様だ。

• sox-14.4.1をダウンロードして、インストール

% gunzip -c sox-14.4.1.tar.gz | tar -xvf -
% cd sox-14.4.1
% ./configure --prefix=/home/username
% make
% make install

• soxのオプション

• wav→mp3変換

# sox dds1000x.wav dds1000x.mp3

• ファイル形式は、拡張子で指定出来る（typeで指定も可能）
• .datはテキスト形式みたいだが、余計なものが付いているので、ちょっと使いづらいかな

音で確認

• 可聴周波数は20kHzぐらいまでなので、サンプリング周波数は1MHzではオーバースペックだ
• 例えばCDだと44.1kHzなので、キリの良いところでサンプリング周波数は50kHzとして計算してみる

DDSの音

• 改めて、DDSによる1kHzの正弦波の生成

% dds 1000 50000 | add 1 | mul 128 | int > dds1000x.dat
% txt2wav -50000 dds1000x.dat
% lame dds1000x.wav dds1000x.mp3
% dds 1000 50000 | add 1 | mul 128 | int | head -250 | gp_wav > dds1000x.png
% dds 1000 50000 | add 1 | mul 128 | int | fftw | gp_fft > dds1000x_fft.png

• 8bitなので、少し高調波もあるが、十分キレイな音に聞こえる
  • ちょっと、こもったような音に感じるかな

三角波の音

• 1kHzの三角波の生成

% triangle 1000 50000 | add 1 | mul 128 | int > tri1000x.dat
% txt2wav -50000 tri1000x.dat
% lame tri1000x.wav tri1000x.mp3
% triangle 1000 50000 | add 1 | mul 128 | int | head -250 | gp_wav > tri1000x.png
% triangle 1000 50000 | add 1 | mul 128 | int | fftw | gp_fft > tri1000x_fft.png

• 正弦波に近い音だが、聞き比べると高調波が多いことが分かる
  • こもった感が少ないので、正弦波よりはっきりした感じ

三角波の音（改）

• 三角波に、偶数次の高調波が含まれているのが気になるなぁ
• あれ、切捨てしなければ出ない。。整数化のための切捨てで非対称になるってことか。。。

  % triangle 1000 50000 | add 1 | mul 128 | fftw | gp_fft > tri1000xx_fft.png
  

• じゃ、こうすれば良いのか！
• 25kHzに、ちょっとゴミ出てるけど、切り捨てない場合と大差無い感じだ

  triangle 1000 50000 | mul 128 | int | add 128 | fftw | gp_fft > tri1000xxx_fft.png
  

• つまり1kHzの三角波の生成は、この方が良さそう

% triangle 1000 50000 | mul 128 | int | add 128> tri1000xxx.dat
% txt2wav -50000 tri1000xxx.dat
% lame tri1000xxx.wav tri1000xxx.mp3
% triangle 1000 50000 | mul 128 | int | add 128 | head -250 | gp_wav > tri1000xxx.png
% triangle 1000 50000 | mul 128 | int | add 128 | fftw | gp_fft > tri1000xxx_fft.png

• 音的には、、あまり変わんないね。。。

矩形波の音

• 1kHzの矩形波の生成
• 2x128=256なので、255を超えないように、127を掛けている

  % square 1000 50000 | add 1 | mul 127 | int > squ1000x.dat
  % txt2wav -50000 squ1000x.dat
  % lame squ1000x.wav squ1000x.mp3
  % square 1000 50000 | add 1 | mul 127 | int | head -250 | gp_wav > squ1000x.png
  % square 1000 50000 | add 1 | mul 127 | int | fftw | gp_fft > squ1000x_fft.png
  

• 矩形波は、かなり高調波が強く感じられる
  • 非常にはっきりとした音という印象

ノコギリ波の音

• 1kHzのノコギリ波の生成

% saw 1000 50000 | add 1 | mul 128 | int > saw1000x.dat
% txt2wav -50000 saw1000x.dat
% lame saw1000x.wav saw1000x.mp3
% saw 1000 50000 | add 1 | mul 128 | int | head -250 | gp_wav > saw1000x.png
% saw 1000 50000 | add 1 | mul 128 | int | fftw | gp_fft > saw1000x_fft.png

• 奇数倍だけでなく偶数倍の高調波が含まれるので、
  • いろいろな音が混ざったような印象がある

ノイズの音

• ノイズの生成

% rand | head -50000 | add 1 | mul 128 | int > randx.dat
% txt2wav -50000 randx.dat
% lame randx.wav randx.mp3
% rand | head -50000 | add 1 | mul 128 | int | head -100 | gp_wav > randx.png
% rand | head -50000 | add 1 | mul 128 | int | fftw | gp_fft > randx_fft.png

• ホワイトノイズっぽい感じになっている

PN符号の音

• PN符号の生成

% pngen | bit0 | head -50000 | mul 255 > pngen0x.dat
% txt2wav -50000 pngen0x.dat
% lame pngen0x.wav pngen0x.mp3
% pngen | bit0 | head -50000 | mul 255 | head -100 | gp_wav > pngen0x.png
% pngen | bit0 | head -50000 | mul 255 | fftw | gp_fft > pngen0x_fft.png

• うぁっ、これは、ちょっと機械のノイズっぽい感じ。全然違うなぁ

DDSの音（サンプリング1MHz）

• WAVファイルを作るために、サンプリング周波数を50kHzにしていたが、
• 考えてみると、サンプリング周波数1MHzのままで、変換できるかもしれない。

• ということで、1MHzサンプリングでWAVファイルを作って、
• lameでMP3に変換すると、ちゃんと48kHzにリサンプリングされた！

  % dds 1000 | add 1 | mul 127 | int > dds1000xx.dat
  % txt2wav -1000000 dds1000xx.dat
  % lame dds1000xx.wav dds1000xx.mp3
  % dds 1000 | add 1 | mul 127 | int | head -5000 | gp_wav > dds1000xx.png
  % dds 1000 | add 1 | mul 127 | int | fftw | gp_fft > dds1000xx_fft.png
  

• 元のサンプリング周波数が大きくても、高調波はあまり改善されているようには見えない。

• 三角波と同じように処理してみた
• この方が上下の対称性ということでは良い筈だけど。。

  % dds 1000 | mul 127 | int | add 128 > dds1000xxx.dat
  % txt2wav -1000000 dds1000xxx.dat
  % lame dds1000xxx.wav dds1000xxx.mp3
  % dds 1000 | mul 127 | int | add 128 | head -5000 | gp_wav > dds1000xxx.png
  % dds 1000 | mul 127 | int | add 128 | fftw | gp_fft > dds1000xxx_fft.png
  

• 高調波的には、うーん。むしろ悪くなってるかも³
  • 非対称ということより、全体的な歪の影響の方が大きいのかな

• まぁ、元のサンプリング周波数が大きくても、特に音質が良くなるわけでも悪くなるわけでも無い様だ。
• でも、サンプリング周波数に関係なく、WAVファイルは作れるし、MP3変換時にリサンプリングされる

• わざわざ、CDで使われる44.1kHzとか、DVD規格の48kHzに近い、50kHzサンプリングで生成していたが、
• 周波数を全く気にせずに、いきなり変換しても全然問題なさそうだ。

２値化

• 信号の２値化の方法について、いろいろ試してみる

単純な２値化

• 単に、コンパレータを通すことで信号を2値化する（150608追記）
• プログラムで生成した正弦波なので、単純に２値化しても、単に矩形波になるだけだと思っていたんだけど。。

  % dds 1000 | comp 0 | head -2000 | gp_wav > comp1000.png
  

• fftwを掛けてみると
• 単純な矩形波なので、奇数倍の高調波のみ含まれるはずだが。。

  % dds 1000 | comp 0 | fftw | gp_fft > comp1000_fft.png
  

• 実際には、偶数倍の高調波もかなり含まれている
  • これは、2値化の際にDutyがずれるのが理由かもしれないが、全部同じ大きさなのは何故だろう？

• 考えてみると、全周波数に均等に分布しているということから、これは量子化誤差だろう。
• 2値化対象の正弦波が計算値なので、キレイに均等になっているのが分かるということか。。

• 量子化誤差が、これほどハッキリと見えるものだとは思わなかった。
• でも逆に言えば、他のどのような２値化の方法でも、この量子化誤差の影響はあるということだろう。

• 音として聞いてみる

% dds 1000 | comp 0 | mul 255 > comp1000.dat
% txt2wav -1000000 comp1000.dat
% lame comp1000.wav comp1000.mp3

• 音としては、普通の矩形波と区別はつかないなぁ。。量子化誤差はレベルが小さいので、ほとんど影響していない。
  • もちろん、元々の正弦波の音ではない。

PWM

• 三角波と比較することによって、信号をPWM化する

• 1kHzの正弦波を、10kHzの三角波と比較することで、PWM化してみた

  % triangle 10000 > tri10000.dat
  % dds 1000 | comp tri10000.dat | head -2000 | gp_wav > pwm1000.png
  

• fftwを掛けてみる

% dds 1000 | comp tri10000.dat | fftw | gp_fft > pwm1000_fft.png

• 肝心の1kHzと同じぐらいのレベルの雑音が含まれるので、このままではマズイ
• 10kHzの三角波でPWM化しているので、その半分の5kHz以下を使えば良さそうだ

• PWM波形を、音で聞いてみよう

% dds 1000 | comp tri10000.dat | mul 255 > pwm1000.dat
% txt2wav -1000000 pwm1000.dat
% lame pwm1000.wav pwm1000.mp3

• FFTの結果を見ると、ローパスフィルタを通して使う必要がありそうだけど、
• ２値なのに、意外とちゃんと聞ける。まぁ、ノイズは結構高い周波数だからかな
• でも、音質はかなり変わっちゃってる感じだ

PWMの周波数を上げる（161223追記）

• フィルタレスのD級アンプの音を聞いて、PWMの周波数を上げれば音質は良くなりそうなので試す
• まず、40kHzの三角波と比較するようにして、fftwを掛けてみた。

  % triangle 40000 > tri40000.dat
  % dds 1000 | comp tri40000.dat | fftw | gp_fft > pwm1000x_fft.png
  

• 雑音の周波数は少し上がっている。大きなのは20kHz以上だし、ちょっと良くなったかも。

• では早速、音を聞いてみよう

% dds 1000 | comp tri40000.dat | mul 255 > pwm1000x.dat
% txt2wav -1000000 pwm1000x.dat
% lame pwm1000x.wav pwm1000x.mp3

• fftの結果では良さそうだけど、音質はあまり良くないなぁ。。
• うーん、じゃ実際のD級アンプと同じぐらいの300kHzの三角波と比較してfftwを掛けてみよう

  % triangle 300000 > tri300000.dat
  % dds 1000 | comp tri300000.dat | fftw | gp_fft > pwm1000xx_fft.png
  

• 大きな雑音は、高い周波数になったけど、全体的に見るとむしろ悪くなった⁴ような。。

• あまり期待できそうもないが、一応音を聞いてみよう

% dds 1000 | comp tri300000.dat | mul 255 > pwm1000xx.dat
% txt2wav -1000000 pwm1000xx.dat
% lame pwm1000xx.wav pwm1000xx.mp3

• あれ、ノイズがほとんど聞こえなくなった。。。
• fftで見ると雑音だらけなんだけど、人間の耳には聞こえないのか。。
• スピーカー自体が、ローバスフィルタの役目をしているのかもしれないなぁ。

• でも、20kHz以上の雑音はあまり影響ないと思っていたんだけど、やっぱり影響するのかな
• まぁ、いずれにしてもPWMの周波数は、けっこう高めじゃないとダメらしい。
• PWMの周波数が上がると、三角波自体の品質というか解像度も影響ありそう

RPWM

• 同様に、乱数と比較することで、信号をPWM化する
• やはり1kHzの正弦波を、乱数と比較してPWM化してみた
  • rand出力は0〜1なので、２倍して１引くことで、-1〜1の範囲に変更している

    % rand | mul 2 | add -1 > rand.dat
    % dds 1000 | comp rand.dat | head -2000 | gp_wav > rpwm1000.png
    

• fftwを掛けてみる

% dds 1000 | comp rand.dat | fftw | gp_fft > rpwm1000_fft.png

• ノイズだらけなんだけど、肝心の1kHzの信号以外のレベルは小さく抑えられている
• 但し、ノイズは全ての周波数に渡っているので、フィルタとかで削るのは難しい
• でも意外と悪くない。そのままでもそれなりに使えるかもしれない。

• RPWM波形も、音で聞いてみよう

% dds 1000 | comp rand.dat | mul 255 > rpwm1000.dat
% txt2wav -1000000 rpwm1000.dat
% lame rpwm1000.wav rpwm1000.mp3

• うーん、音質自体はそんなに悪くないけど、ホワイトノイズは、そのまま残ってるなぁ。。
  • やっぱり聞いてみると、問題点も分かりやすいね。

１次ΔΣ

• １次のデルタシグマによる２値化

• これをそのままプログラムにすると、こんな感じかな（ds1.c）
• 面倒なので、入力は-1〜1を想定、積分器はdoubleにしました

  #include <stdio.h>
  
  int     main(int argc, char *argv[])
  {
    int    out;
    double in, delta, sigma=0;
  
    while (scanf("%lf", &in) != EOF) {
      delta = in - out;
      sigma += delta;
      out = (sigma > 0) ? 1 : -1;
      printf("%d\n",out);
    }
    return 0;
  }
  

• 先頭の2000ポイントのみプロットしてみる

% dds 1000 | ds1 | head -2000 | gp_wav > ds1_1000.png

• FFTを掛ける

% dds 1000 | ds1 | fftw | gp_fft > ds1_1000_fft.png

• いわゆるノイズシェーピングにより、ノイズ成分が高周波に押し上げられていることが分かる

• 音を聞いてみる

% dds 1000 | ds1 | add 1 | mul 127 | int > ds1_1000.dat
% txt2wav -1000000 ds1_1000.dat
% lame ds1_1000.wav ds1_1000.mp3

• ノイズシェーピングの効果は絶大！ ほとんどのノイズが可聴周波数外になっているらしい
  • とても２値の信号とは思えないほど、ノイズが少なく感じる

• ひょっとすると、MP3へのリサンプリング時にLPFが掛かっている可能性もあるけど

擬似ホワイトノイズを作る

• 前述の項目で、rand関数を使ったホワイトノイズを作ったんだけど、
• やはり乱数だけに、周波数毎の大きさに結構バラツキがある。

• 特に、デジタルフィルタの特性を見たい時には、バラツキが無いほうが良いなぁ。
• ということで、fftwの実験用の擬似ホワイトノイズを作ってみる

正弦波の加算プログラム

• 理屈で言えば、全ての周波数の正弦波を足していけば、ホワイトノイズになる
• ddsで正弦波を生成して、addで加算すれば良いが、とても面倒なのでプログラムを作る

• DDSで、1Hzおきに正弦波を生成したい。
• DDSで生成可能な周波数は、サンプリング周波数/波形テーブルの大きさの整数倍になるので、
• 波形テーブルは、少なくともサンプリング周波数分だけ必要になる

• 最初は、試しにサンプリング周波数1kHzで、0〜999Hzを加算してみた
• すると、結果が常に0になってしまった。プログラムを何度見直してもゼロになる。

• よく考えてみると、ナイキスト周波数で折り返しになっているので、全部足せば０なのか。。
• ということで、サンプリング周波数の半分まで加算することにした。（white.c）

  #include<stdio.h>
  #include<stdlib.h>
  #include<math.h>
  
  #define MAX_SAMPLE (1000000)
  #define SAMPLE (1000)
  
  double sin_table[MAX_SAMPLE/4+1];
  int t[MAX_SAMPLE];
  int N, f1, f2;
  
  void make_table(void)
  {
          int x;
  
          for (x=0; x<=N/4 ;x++)
            sin_table[x] = sin(2*M_PI/N*x);
  }
  
  double sin_x(int x)
  {
          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, j, d;
          double out;
  
          N = SAMPLE;
          if (argc>1)
            N = atoi(argv[1]);
          f1 = 0;
          f2 = N/2-1;
          if (argc>2)
            f2 = atoi(argv[2]);
          if (argc>3)
            f1 = atoi(argv[3]);
          make_table();
  
          srand(777);
          for (i=0; i<N ;i++) t[i] = rand()%N;    // 位相の初期値を０ではなく、バラつかせる
  
          for (i=0; i<N ;i++) {
            out = 0;
            for (j=f1; j<=f2 ;j++) {
                  out += sin_x(t[j]);
                  d = j;
                  t[j] = (t[j]+d)%N;
            }
            printf("%f\n",out);
          }
  }
  

• サンプリング周波数 Sと、生成する周波数範囲 f1〜f2 の指定も可能
  • パラメータは、S, f2, f1の順に指定する
  • 周波数範囲を省略すると、1Hz間隔で 0〜S/2の範囲で生成する

サンプリング周波数1kHz（失敗）

• サンプリング周波数１kHzで、加算した結果をFFTで見てみる

% white0 1000 | fftw | gp_fft > white01000_fft.png

• よしよし、思った以上に、ちゃんと加算されてるなぁ。

• では、波形を見てみよう
• 先頭の100ポイントのみプロットしてみる

  % white0 1000 | head -100 | gp_wav > white01000_100.png
  

• えぇー、なんだこれ。。

• 波形の全体をみてみる

% white0 1000 | gp_wav > white01000.png

• うーん。fftwの実験用には、これでも良いんだけど。。
  • ホワイトノイズには程遠いなぁ。。それにちょっと使いづらい。。。

• これは、ひょっとすると δ関数⁵を作ってるようなものなのかもしれない。。

位相特性用には使える

• 上で作成した変な波形も、位相が揃っているので、デジタルフィルタの位相特性を見るためには使える（140905追記）
• この信号の位相特性をみてみると⁶、

  % white0 1000 | fftw | fftw_unwrap | gp_ang > white01000_ang.png
  

• ナイキスト周波数まで、一定の位相（π、或いは-π）になっている。
  • でも、ナイキスト周波数以上は、含まれないので、位相はいきなり０となる（不連続になっている）
  • これは、FFTではそういうものなのか、fftwの仕様なのか分からないが、ナイキスト周波数以上は使えない様だ

• この信号を、例えばRCフィルタを通してみると、

% white0 1000 | rcf 0.5 | fftw |fftw_unwrap | gp_ang > rcf_ang.png

• RCフィルタの位相特性が見える（でも±が逆のような）⁷

• まぁ、いずれにしても何かの役に立つかもしれないので、
• 擬似ホワイトノイズ０として、使えるようにしておく（white0.c）

サンプリング周波数1kHz（改）

• 多分、正弦波を合成する時に、スタート位置が一緒なので、
• 全部０から始まって、ピークの位置が一致してしまうんだろう

          for (i=0; i<N ;i++) t[i] = 0;
  

• ということで、各周波数の位相の初期値を０じゃなくて、
• 乱数でバラつかせることにした。

          srand(777);
          for (i=0; i<N ;i++) t[i] = rand()%N;
  

• では、改めて先頭の100ポイントのプロットと、FFT結果を見てみよう

% white 1000 | head -100 | gp_wav > white1000_100.png
% white 1000 | fftw | gp_fft > white1000_fft.png

• 今度は、大丈夫そうだ。

• 初期値０の場合と、最小値・最大値を比較してみる

% white0 1000 | minmax
-318.308839, 318.308839

% white 1000 | minmax
-51.121728, 45.900439

• 波形の幅もかなり抑えられていることがわかる。
  • ピークの位置がうまく分散されているようだ。

サンプリング周波数48kHz

• では、音で聞いてみることを前提に、サンプリング周波数を48kHzにして
• 正弦波を加算してみる

• やはり先頭の100ポイントのプロットと、FFT結果

% white 48000 | head -100 | gp_wav > white48000_100.png
% white 48000 | fftw | gp_fft > white48000_fft.png

• 素晴らしい！少なくともfftwでみる限りは、全周波数について同レベルの信号になっている
  • あくまで計算上なので、実際の波形にする場合は、そんなにうまく行くわけない⁸と思うけど

• 最小値・最大値を確認する

% white 48000 | minmax
-451.917214, 435.368105

• 16bitのWAVファイルは、-32768〜32767のデータ範囲なので、
• 32768/452≒72.5ぐらいだから、ｘ70ぐらいでいいか。

• WAVファイルを生成、MP3に変換する

% white 48000 | mul 70 > white48000.dat
% txt2wav -48000 white48000.dat
% lame white48000.wav white48000.mp3

• ちゃんと、ホワイトノイズっぽい感じになってる

• まぁ、ほんとにホワイトノイズかどうかは別にしても
• fftwの実験用には、乱数によるノイズよりも、ずっと使いやすいと思う

指定範囲だけ生成

• 指定範囲だけ、生成できる様にしてみた（140902追記）

% white 1000 100 10 | fftw | gp_fft > white1000_100_10_fft.png

• サンプリング周波数1kHzで、10〜100Hzの範囲だけ生成してみた
  • 少しノイズもあるけど、レベルは十分小さいので、全然問題なし⁹

• 乱数だと、こんなことは無理だけど、これは加算してるだけなので簡単だ
• 実験用のデータとしては、結構使い道があるんじゃないかと思う

• 関係ないけど、一応、音にしてみた

% white 1000 100 10 | mul 5 | add 128 > white1000_100_10.dat
% txt2wav -1000 white1000_100_10.dat
% lame white1000_100_10.wav white1000_100_10.mp3

• はっきりしなくて、音とは言えない感じだけど、こもったような感じでちょっと面白い。
• 何かの効果音に使えそうだ

まとめ

• プログラムで波形を作っても、波形を見ただけでは違いが分からなかったりするけど、
• fftwがあれば、なんでもFFTを掛けて、すぐに含まれる周波数を確認できる

• Unix環境なら、gccでプログラムを作って、
• コマンドラインで、フィルタを組み合わせて、いろいろ実験できる。
• プログラムでデジタルフィルタでは、いろいろ試してみた

• サンプル数の制限も無いので、いろいろ試せて、とっても便利
• 個人的には、デジタル系の処理を考えるときに、fftwは無くてはならない道具になっている。

参考

• C言語による乱数生成
• 乱数生成も、いろいろ奥が深そうだ

• WAVファイルフォーマット
• txt2wavプログラムを、使わせてもらいました。感謝！

• ΔΣ変調の部屋
• ΔΣ変調について、詳しく解説されている

• 始めにΔ変調ありき
• ΔΣ変調が考案された経緯が分かって面白かった