音についてまとめてみる（３） - 「がんばれない」けど「がんばりたい」

フーリエ変換をコンピュータで取り扱う時は、離散フーリエ変換からFFTに話が進むのだが、

用語を本やネットで調べれば調べるほど、混乱したので、

用語のまとめ
プログラムによる計算
今回vDSPを使用した結果

についてまとめようと思う。

詳しい事は、がんばってネットなどで調べて頂きたいのですが、自分の理解としては、以下の様になった。

フーリエスペクトル-------

離散フーリエ変換を行うと結果は、実部と虚部に分かれる。この２部を用いた表現。

X(k) = Real(k) + jImg(k)

１．パワースペクトル

振幅スペクトルとも言う（と思う）。実部を横軸、虚部を縦軸とした時、

パワースペクトル＝√Real(k)Real(k) + jImg(k)jImg(k) と計算する。

※ピタゴラスの定理検証プログラムで以下の条件でsin波を生成し、 FFTを行い、パワースペクトルを計算。そしてパワースペクトルからdB値を計算してみた。

基本周波数：2Hzの音を1秒間サンプリング周波数：16Hz

//基本周波数
double freq = 2.0;

//サンプリング周波数
double sampleRate = 16.0;

//各サンプルポイントの角度を決める変数
double phase = 0.0;

//各サンプルポイント間の角度（ラジアン）間隔
double unit = freq * 2 * M_PI / sampleRate;

//音をながす秒数
float sec = 1.0;

//総サンプルポイント数
int sampleLen;
sampleLen = (int)sec * sampleRate;

//
//FFTを行うのでreal, imagの配列を用意
//2Hzの音を作成
//
float real[sampleLen], imag[sampleLen];
for (int i = 0; i &lt; sampleLen; i++) {
real[i] = sin(phase);
phase += unit;
// 虚部は0
imag[i] = 0.0;
}

//FFT
//MyFFTクラスではハニング窓を使用し、vDSP_fft_zripでFFTをしています
MyFFT* myfft = [[MyFFT alloc] initWithCapacity:sampleLen];
[myfft process:real];

//
//パワースペクトル
// powerSpectrum配列にパワースペクトル値が格納される
float powerSpectrum[sampleLen/2];
vDSP_vdist(myfft.realp, 1, myfft.imagp, 1, powerSpectrum, 1, sampleLen/2);

//パワースペクトルはサンプル数で割る
for (int i = 0; i &lt; sampleLen/2; i++) {
powerSpectrum[i] = powerSpectrum[i] / (sampleLen/2);
}

//
// dB計算
//
float zeroReference[sampleLen/2];
float decibelChannelL[sampleLen/2];

// 1.0 = 0dBとする
for (int i = 0;i&lt;sampleLen/2;i++)
{
zeroReference[i] = 1.0;
}

// vDSPでのdB計算
vDSP_vdbcon (
powerSpectrum,　//入力
1,
zeroReference, //ゼロリファレンス
decibelChannelL,//結果を格納する配列
1,
sampleLen / 2, //サンプル数
1 //0:入力をパワーとして計算。1:amplitudeとして計算。今回はamplitude（振幅）
);

//結果を表示
for (int i = 0; i &lt; sampleLen/2; i++) {
float s = powerSpectrum[i];
float db = 20*log10(s);
printf("%dHz %f %f %fn", i, s, db, decibelChannelL[i]);
}

これで、20log10(power)とvDSPでの計算結果が同じになった。