cで書き出したwavファイルの周波数が思う結果と違う。

Question

C言語で音声合成もどき ～WAVファイルを生成する～

を参考に、指定した時間（秒）の内に指定した開始時の周波数から終了時の周波数まで音の高さが変わるコードを書こうとしています。

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

// 44100Hz, 8bit, 1ch のWAVデータ
#define SMPL 44100
#define BIT  8

void wav_write(const char *filename, char *buf, int size)
{
    int filesize = 44 + size;
    char *work = (char *) malloc(filesize);
    FILE *fp = fopen(filename, "wb");

    if (fp == NULL) return;

    /* RIFFヘッダ */
    memcpy(work, "RIFF", 4);
    work[4] = (filesize - 8) >> 0  & 0xff;
    work[5] = (filesize - 8) >> 8  & 0xff;
    work[6] = (filesize - 8) >> 16 & 0xff;
    work[7] = (filesize - 8) >> 24 & 0xff;
    /* WAVEヘッダ */
    memcpy(work+8, "WAVE", 4);
    /* fmtチャンク */
    memcpy(work+12, "fmt ", 4);
    work[16] = 16;
    work[17] = work[18] = work[19] = 0;
    work[20] = 1;
    work[21] = 0;
    work[22] = 1;
    work[23] = 0;
    work[24] = SMPL >> 0  & 0xff;
    work[25] = SMPL >> 8  & 0xff;
    work[26] = SMPL >> 16 & 0xff;
    work[27] = SMPL >> 24 & 0xff;
    work[28] = (SMPL * (BIT / 8)) >> 0  & 0xff;
    work[29] = (SMPL * (BIT / 8)) >> 8  & 0xff;
    work[30] = (SMPL * (BIT / 8)) >> 16 & 0xff;
    work[31] = (SMPL * (BIT / 8)) >> 24 & 0xff;
    work[32] = ((BIT / 8)) >> 0 & 0xff;
    work[33] = ((BIT / 8)) >> 8 & 0xff;
    work[34] = BIT >> 0 & 0xff;
    work[35] = BIT >> 8 & 0xff;
    /* dataチャンク */
    memcpy(work+36, "data", 4);
    work[40] = size >> 0  & 0xff;
    work[41] = size >> 8  & 0xff;
    work[42] = size >> 16 & 0xff;
    work[43] = size >> 24 & 0xff;
    memcpy(work + 44, buf, size);

    /* 書き出し */
    fwrite(work, filesize, 1, fp);
    fclose(fp);
    free(work);
}

int main(void)
{
    int size;
    char *buf;
    int sec;

    sec = 3;
    size = (SMPL*(BIT/8)) * sec;
    buf = (char *) malloc(size);

    double freq_start = 886;        // 開始時の周波数
    double freq_end = 443;         // 終了時の周波数

    double t;        // 時間（秒）
    double freq;    // 周波数
    double y;        // 出力

    for (int i = 0; i < size; i++) {
        // 時間（秒）
        t = i / (double)size * sec;

        // 周波数
        freq = freq_start + (t / sec) * (freq_end - freq_start);

        // 出力
        y = sin(2.0 * M_PI * freq * t);

        buf[i] = floor(y * 63) + 128;

#ifdef EBUG
        printf("%f\n", freq);
#endif
    }

    wav_write("test.wav", buf, size);

    free(buf);
    return 0;
}

この例では周波数を886Hzから443Hzまで変化させようとしていますが、
実際に出力されたtest.wavは886Hzから0Hzまで落ちてしまいます。

ちなみに

    double freq_start = 886;        // 開始時の周波数
    double freq_end = 443;         // 終了時の周波数

を

    double freq_start = 443;        // 開始時の周波数
    double freq_end = 886;         // 終了時の周波数

とすると意図した通り周波数は443Hzから886Hzに変化します。

何がいけないのでしょうか？

(補足)

printfで出力すると周波数(freq)は思い通りに変化しています。
（自分だけじゃなく、Zuishinさんがダブルチェックをしていただけたのでこれは大丈夫）

所が音源ファイル(test.wav)を実際に聴いてみると、
周波数は0まで落ちているのです。
その確認は「audacity」という、波形を見ながら音を再生できるソフトで出来ますし、
耳で聞いてもおかしいのは分かります。

周波数は正しく変化していて、
その上でwavファイルがその通りに書き出されてない。

freqとtest.wavが一致しないのは何故なのでしょう？

Answer 1

原因は位相の求め方にあります。

本件ではサンプルごとになだらかに周波数を変えようとしていますね。このときある瞬間の周波数は直前の位相との位相差により決まります。ところが質問コードでは隣り合ったサンプルの位相の差は期待した周波数になりません。そのため意図した周波数とずれた周波数の音になってしまうのです。

y = sin(2 * M_PI * freq * t);

で波形を求めるためには条件が必要で、その条件はこの式を用いる期間で「周波数が不変でt=0の時の位相が2πの倍数」であることです。

以下のようにするとサンプル間の位相差がi番目のサンプルの周波数に従った値になるため期待通りの結果になります。

...
double phase = 0.0; // 位相
for (int i = 0; i < num_sample; i++) {
  double freq = ...;
  phase += 2 * M_PI * freq / SMPL;
  double y = sin(phase);
  ...
}
...

---

他の指摘：

• forループの繰り返し回数

本来の意味は「サンプル数」であるべきです。しかし、コードでは「バッファーのバイト数」となっています。本コードではたまたま両者が同じ値なので問題は生じませんが正しいコードとは言えず、例えばこの論理を16bit PCMへ流用した場合などにバグの元になると思います。

Comments

[taro_nii_chan]

なるほど！
って言えたらかっこいいんでしょうけど、
理論がまだ良くわかってません。
この辺のことはどこで勉強できますか？
書籍でもいいしサイトでもいいので
参考になる物があったら教えて下さい。

取り敢えず、思い通りじゃなかったものが動いて
すごく爽快な気分です。

ありがとうございます。

[KSwordOfHaste]

人それぞれ気づき方は違うかも知れません。自分は円周上を角速度をだんだん遅くしながら回転する運動を頭に思い描き「直前の点の位置(角度)からどのくらい角度が変化するかが周波数になるのだから位相が連続してないといけないんだな・・・」というふうに気づけました。

勉強という意味では音波＝波についての数学（高校や大学初等の数学）や微分的な考え方をよく理解しておくことが役立ちます。理解というのは数式や公式を単に覚えるのではなく幾何的・物理的な意味のイメージを持てるようにすることです。専門家しかわからないような難しい話だと自習は困難ですが、高校などで学んだはずの数学が基礎になっている応用問題ならば、（学校で学んだ当時は公式丸暗記だったとしても）プログラミングで実際に使うとわかったなら「より真剣に深く理解しよう」という気が起きると思います。そのやる気を大事にして関連しそうなサイトを眺めていろいろ考える訓練をすれば、忘れていた数学的思考がよみがえってくると思います。

[taro_nii_chan]

分かりました。
ありがとうございました。
またよろしくお願いします。

[ozwk]

> 質問コードでは隣り合ったサンプルの位相に連続性がありません。

位相が微分可能なので、連続性はあるのでは?
(そもそも時間が離散ですが、離散化前の式が。)

[KSwordOfHaste]

>ozwkさん
ご指摘のように「連続性」という言葉は不適切でした。

i番目のサンプルの周波数がfならi-1番目のサンプルの位相p-Δとi番目のサンプルの位相pが「周波数fとしての位相の差になってなければならない」というべきでしたね。

それでもtの大きさが大きくなるにつれて位相差のかい離は大きくなっていき連続でないような外観を呈するような気がしたので「連続でない」といっちゃいました。

[KSwordOfHaste]

遅くなりましたが、連続性という用語を回答で使わないよう表現を変更しました。

Answer 2

bufを書き出してみました(図は最後の方だけ)

これ見ればわかるように最後の方で周波数が0になります。

ということは、そもそも根本的に、
質問文にあるようなアルゴリズムで886から443Hzに変化する音は作れません。

というのも、音の瞬間的な周波数(瞬時周波数)は音の波形が
sin(2π*φ(t))で表せるとき、
dφ/dtで表せます。

単純に、sin(2π*f(t)*t)だからf(t)が今の周波数になるわけじゃないです。

これに基づいて計算すると、
今回の波形の瞬時周波数は
fs + 2*(fe-fs)*t/Tとなり、(t=0~T)
t=Tのとき
2fe-f0 = 2443 - 883 =0となります。

逆に、意図通りの周波数変化:f(t)となる音を出すには、
f(t)の積分F(t)を計算して、
sin(2π*F(t))とすれば良いことがわかります

Comments

[taro_nii_chan]

ありがとうございます。
質問に追記しましたが、
「audacity」というソフトでこの様な波形が見られます。
実際、ozwkさんが貼っていただけた画像と同様のものが現れています。

[ozwk]

なので、y=...の式がそもそもおかしいんです、と言いたかったんですが
追記します。

[taro_nii_chan]

素晴らしい。納得です。
KSwordOfHasteさんのアプローチも素晴らしくて
ベストアンサーを1人しか決められないのが恨めしいくらいです。

Answer 3

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

int main(void) {
    int size;
    int sec;

    sec = 3;
    size = 30;

    double freq_start = 886;        // 開始時の周波数
    double freq_end = 443;         // 終了時の周波数

    double t;        // 時間（秒）
    double freq;    // 周波数

    for (int i = 0; i < size; i++) {
        // 時間（秒）
        t = i / (double)size * sec;

        // 周波数
        freq = freq_start + (t / sec) * (freq_end - freq_start);

        printf("%f\n", freq);
    }
    return 0;
}

上記の出力結果は次のようになりました。

886.000000
871.233333
856.466667
841.700000
826.933333
812.166667
797.400000
782.633333
767.866667
753.100000
738.333333
723.566667
708.800000
694.033333
679.266667
664.500000
649.733333
634.966667
620.200000
605.433333
590.666667
575.900000
561.133333
546.366667
531.600000
516.833333
502.066667
487.300000
472.533333
457.766667

i < size という条件なので 443 ちょうどにはなりませんが、886 から 443 へ変化しています。

Comments

[taro_nii_chan]

そうなんです、printfで出力すると周波数(freq)は思い通りに変化しているのです。

所が音源ファイル(test.wav)を実際に聴いてみると、
周波数は0まで落ちているのです。
その確認は「audacity」という、波形を見ながら音を再生できるソフトで出来ますし、
耳で聞いてもおかしいのは分かります。

周波数は正しく変化していて、
その上でwavファイルがその通りに書き出されてない。
全く不思議でしょうがないのですが。

これは何なんでしょう？

[Zuishin]

その情報は質問へ追記してください。freq が正しいのであれば、freq の値を追う作業は無意味でした。

[taro_nii_chan]

質問への追記、しました。

freqは「私の確認した限りでは」だったので、
そもそもそこからしておかしいという事も考えられたので、
freqを追う作業は私にとっては全然無意味ではないです。
ありがとうございます。