Cでフーリエ変換してパワースペクトル密度を求めるプログラムを作成する

C

ちょっと諸事情により、C言語で音声データのパワースペクトル密度を求めないといけなくなってしまったのでそのメモ書き（どんな事情

パワースペクトル密度とは

まず、パワースペクトル密度とはなんでしょうか（全くわかりません）

一応Wikiのパワースペクトル密度のページも読みましたが、何を言ってるんだこいつは状態でした。

これについて知るには、まずパワースペクトルがなんなのかについて一度まとめておきます。

http://www.tagen.tohoku.ac.jp/labo/ishijima/FFT-02.html

フーリエ変換もパワースクペクトルも、グラフに表すと横軸は周波数になるみたいです。そして縦軸の値が異なる。

フーリエ変換は縦軸を波形の周波数成分を振幅として表したもの、パワースペクトルは縦軸を波形の周波数成分をエネルギーとして表したものとのことです。

そして単位は以下のようになります。

フーリエ変換：[m] パワースペクトル：エネルギー：[m2/Hz]，もしくはその平方根

そしてパワースペクトル密度とはパワースペクトルとほぼ同義で使われているようです。

強いて言うなら、言葉の意味的に以下のような違いな気がします。

うーん、なんか、曖昧ですが...笑

まぁやってみないことにはわからない、とりあえずプログラムを作って走らせてみます。

Cでパワースペクトル密度を求める

パワースペクトル密度を求めるプログラム

音声データからパワースペクトル密度P(f_k)をC言語のプログラムにより計算し、そのグラフを描きなさい。ただしk=0, 1, ..., N/2-1 に対し、横軸にf_k = k/NΔ [Hz], 縦軸に10log_10_P(f_k) [dB] をプロットする。

（離散）フーリエ変換については以下のブログがすごくよく書かれてて楽しい。ただしC++なのでauto&やconst auto&は使えないか。

プログラムの国のフーリエ変換 - nursの日記

  • プログラムの世界でできるのは離散フーリエ変換のみ (微積分みたいな近似的なことができないから？）

また以下のプログラムも直感的でわかりやすいと思う。

フーリエ変換とスペクトルの強さの数値計算例(C言語)

そして最終的に高速フーリエ変換でやりたいと考えてたので、同じブログの以下の記事も参考に。

高速フーリエ変換(FFT)をおじさんもC++で作ってみたよ - nursの日記

うむ、高速フーリエ変換の式の導出に関しては1mmも理解できるとは思えないけど...笑

プログラムとしてもすごく長いみたいだけど、sinとcosがなくても計算できる+再帰だから速いのかな。

ソースコードは以下のようになった。フーリエ変換で求めた数値を、spectram.datというファイルに書き込む設定になっている。

power-spectral-density-calculator/dft.c at master · totzyuta/power-spectral-density-calculator · GitHub

また、今回の形式の音声データをパースして正規化するスクリプトもRubyで実装した。

power-spectral-density-calculator/normalize.rb at master · totzyuta/power-spectral-density-calculator · GitHub

グラフに描写

これをグラフに描写していく

gnuplotのインストール

まず、以下を参考にMacでgnuplotを使えるようにしておく。

gnuplotをHomebrewからインストールするときの手順 - Qiita

AquaTermが必要なのでダウンロードして、インストールしておく。　

AquaTerm (Mac OS X graphics terminal) download | SourceForge.net

さきほどダウンロードしたaquatermを指定してbrewでgnuplotをインストール。

$ brew install gnuplot --with-aquaterm

グラフを描写

以下でgnuplotを対話的に操作できる。

$ gnulot

まず最初のグラフは

gnuplot> set xlabel '[Hz]'
gnuplot> set ylabel '[dB]'
gnuplot> plot 'spectrum.dat' with lines

のようにすると、spectrum1.pdfのようなグラフを得ることができた。

また、

gnuplot> set xlabel '[s]'
gnuplot> set ylabel '[V]'
gnuplot> plot 'voice-data.csv' with lines

ようにして自分の音声データをグラフとして描写したところ、voice-graph.pdfとなった。

また、最初に作成したspectrum.datよりピーク時の周波数を100として、２番目に作成した自分の音声データと、周波数100の正弦波のグラフを重ねて表示させてみる。コマンドは以下のように行う。

gnuplot> set xlabel '[Hz]'
gnuplot> set ylabel '[dB]'
gnuplot> plot "data-voice.csv" with lines, 0.01*sin(2*3.14*100*x-8.5) with lines

結果はdata-voice.pdfのようになった。

ピーク時の周波数を自分の音声の測定データから読み取ることはむずかしく、目測で100Hzとしたところ大きな波の動きは似た傾向を示しているように見えるグラフとなった。