WAVファイル解析C（libsndfile、fftw3）

Question

私はWAVファイル内の特定のタイムスタンプで特定の周波数範囲で0〜100の値を与えることができる単純なCアプリケーションを開発しようとしています。

例：44.1kHz（標準MP3ファイル）の周波数範囲があり、その範囲をn個の範囲（0から始まる）に分割したいとします。私は今WAV-のデータを読み取ることができるよ

をLibsndfileが使用：私は、私がこれまで管理してきた

0から100までであること、各範囲の振幅を取得する必要がありますファイル。

しかし、私のFFTの理解はかなり限定されています。しかし、私が必要とする範囲で振幅を得ることが必要であることはわかっています。しかし、私はここからどうやって進むのですか？その目的に合ったライブラリFFTW-3が見つかりました。

私はここにいくつかの助けが見つかりました：https://stackoverflow.com/a/4371627/1141483

を、ここでFFTWのチュートリアルを見て：http://www.fftw.org/fftw2_doc/fftw_2.html

しかし、私はFFTWの行動が不明だとして、私はここから進行するのか分かりません。

もう1つ質問：libsndfileを使用すると仮定します。読み込みを強制的にシングルチャンネル（ステレオファイル）にしてサンプルを読み込むとします。あなたは実際には全体のファイルのサンプルの半分しか読み取っていませんか？それらの半分がチャンネル1から来ているか、自動的にそれらをフィルタリングしますか？

ご協力いただきありがとうございます。

EDIT：私のコードはここで見ることができます：

double blackman_harris(int n, int N){ 
double a0, a1, a2, a3, seg1, seg2, seg3, w_n; 
a0 = 0.35875; 
a1 = 0.48829; 
a2 = 0.14128; 
a3 = 0.01168; 

seg1 = a1 * (double) cos(((double) 2 * (double) M_PI * (double) n)/((double) N - (double) 1)); 
seg2 = a2 * (double) cos(((double) 4 * (double) M_PI * (double) n)/((double) N - (double) 1)); 
seg3 = a3 * (double) cos(((double) 6 * (double) M_PI * (double) n)/((double) N - (double) 1)); 

w_n = a0 - seg1 + seg2 - seg3; 
return w_n; 
} 

int main (int argc, char * argv []) 
{ char  *infilename ; 
SNDFILE  *infile = NULL ; 
FILE  *outfile = NULL ; 
SF_INFO  sfinfo ; 


infile = sf_open(argv [1], SFM_READ, &sfinfo); 

int N = pow(2, 10); 

fftw_complex results[N/2 +1]; 
double samples[N]; 

sf_read_double(infile, samples, 1); 


double normalizer; 
int k; 
for(k = 0; k < N;k++){ 
    if(k == 0){ 

     normalizer = blackman_harris(k, N); 

    } else { 
     normalizer = blackman_harris(k, N); 
    } 

} 

normalizer = normalizer * (double) N/2; 



fftw_plan p = fftw_plan_dft_r2c_1d(N, samples, results, FFTW_ESTIMATE); 

fftw_execute(p); 


int i; 
for(i = 0; i < N/2 +1; i++){ 
    double value = ((double) sqrtf(creal(results[i])*creal(results[i])+cimag(results[i])*cimag(results[i]))/normalizer); 
    printf("%f\n", value); 

} 



sf_close (infile) ; 

return 0 ; 
} /* main */ 

Answer 1

まあそれはすべてあなたが後にしている周波数範囲に依存します。 FFTは、2^n個のサンプルを取り、2 ^（n-1）の実数と虚数を提供します。私は、これらの価値が表現しているものにはかなりぼんやりしていることを認めなければなりません（私は財政上の問題があるときに私が作ったローンの代わりに、円の周りの角度。実際には、オリジナルの2^nサンプルを完全に再構成できる各周波数ビンの正弦と余弦の角度パラメータのアークコスを提供します。

とにかく、これは、実数部と虚数部のユークリッド距離（sqrtf（（real * real）+（imag * imag）））を取ることによって、大きさを計算できるという大きな利点があります。これにより、正規化されていない距離値が得られます。この値を使用して、各周波数帯域の大きさを構築することができます。

したがって、10 FFT（2^10）を注文してください。 1024サンプルを入力します。これらのサンプルをFFTすると、512の虚数値と実数値が戻されます（これらの値の特定の順序は、使用するFFTアルゴリズムによって異なります）。つまり、44.1Khzのオーディオファイルの場合、各ビンは44100/512Hzまたはビンあたり〜86Hzを表します。

これから際立つべきものは、画像などの多次元信号を扱うときに、より多くのサンプル（時間ドメインまたは空間ドメインと呼ばれるもの）を使用すると、より良い周波数表現（周波数ドメイン）。しかし、あなたはもう一方を犠牲にします。これは物事が進む方法であり、あなたはそれで生きなければなりません。

基本的には、必要なデータを得るために、周波数ビンと時間/空間解像度を調整する必要があります。

最初に少しの命名法がありました。以前に言及した1024のタイムドメインのサンプルはあなたのウィンドウと呼ばれます。一般に、この種のプロセスを実行するときは、FFTを使用して次の1024サンプルを取得するために、ある量だけウィンドウをスライドさせたいと思うでしょう。分かりやすいことは、サンプル0から1023、次に1024-> 2047などを取ることです。これは残念なことに最良の結果をもたらしません。理想的には、ウィンドウをある程度重複させて、時間の経過とともによりスムーズな周波数変化を得ることができます。最も一般的には、ウィンドウサイズを半分にスライドします。つまり、最初のウィンドウは0〜> 1023、2番目の512〜> 1535などとなります。

これでさらに問題が発生します。この情報は完全な逆FFT信号再構成を提供しますが、周波数がある程度サラウンドビンに漏れてしまうという問題があります。この問題を解決するために、私よりはるかに知的な数学者がwindow functionのコンセプトを思いついた。窓関数は周波数領域ではるかに優れた周波数分離を提供しますが、時間領域では情報が失われます（つまり、窓関数AFAIKを使用した後に信号を完全に再構築することは不可能です）。

長方形のウィンドウ（事実上何も信号を出さない）から、はるかに優れた周波数分離を提供するさまざまな関数までさまざまなタイプのウィンドウ関数があります（周囲の周波数があなたの興味を引くかもしれません！ ！）。悲しいかな、誰もサイズに合っているわけではありませんが、私はblackmann-harrisウィンドウ関数の大きなファンです（スペクトログラム用）。私はそれが最高の見た目の結果を与えると思います！

しかし先に述べたように、FFTでは非正規化スペクトルが得られます。スペクトルを正規化するには（ユークリッド距離計算後）、すべての値を正規化係数で除算する必要があります（詳細はhere）。

この正規化では、0と1の間の値が得られます。この値を100倍して、0〜100のスケールにすることは簡単です。

ただし、これは終了位置ではありません。あなたがこれから得るスペクトルは、むしろ満足していません。これは、線形スケールを使用してマグニチュードを調べているためです。残念ながら、人間の耳は対数スケールで聞きます。これは、スペクトログラム/スペクトルの見え方にむしろ問題を引き起こします。

これを丸めるには、これらの0から1の値（「x」と呼ぶ）をデシベルに変換する必要があります。標準的な変換は20.0f * log10f(x)です。これにより、1が0に変換され、0が-infinityに変換された値が提供されます。あなたのマグニチュードは適切な対数スケールになりました。しかし、それはいつも役立つとは限りません。

この時点で、元のサンプルビット深度を調べる必要があります。 16ビットサンプリングでは、32767と-32768の間の値が得られます。つまり、dynamic rangeはfabsf（20.0f * log10f（1.0f/65536.0f））または〜96.33dBです。だから私たちはこの価値を持っています。

上記のdB計算から得られた値を取ってください。これに-96.33の値を加えます。明らかに最大振幅（0）は96.33です。今度は同じ値でディディビッドして、今度は、無限大から1.0fまでの値を持っています。下端を0に締めると、0〜1の範囲があり、100を掛けて最終的な0〜100の範囲になります。

これは当初意図していた以上のモンスターポストですが、入力信号に対して優れたスペクトラム/スペクトログラムを生成する方法については十分な根拠があります。

と

さらに（すでにそれを発見したオリジナルポスター以外の人のために）読みを吸わ：

Converting an FFT to a spectogram

編集：さておき、私はキスを見つけたようFFTを使うのがはるかに簡単ですが、フォワードfftを実行するコードは次のようになります：

CFFT::CFFT(unsigned int fftOrder) : 
    BaseFFT(fftOrder) 
{ 
    mFFTSetupFwd = kiss_fftr_alloc(1 << fftOrder, 0, NULL, NULL); 
} 

bool CFFT::ForwardFFT(std::complex<float>* pOut, const float* pIn, unsigned int num) 
{ 
    kiss_fftr(mFFTSetupFwd, pIn, (kiss_fft_cpx*)pOut); 
    return true; 
}