ScilabでFFT(前編) - tera-technology

初記事は大学で出た信号処理の課題についてちょっとやってみたので備忘録ついでに説明しようと思います。

今回のやることは

wavファイルをスペクトル解析して周波数成分を特定
ある周波数をカットするフィルタの作成

です。

早速始めていきましょう

最初にwavファイルのスペクトルを解析します。

こちらがscilabのコードです。

まず解析したいwavファイルを用意します。ここでは sample.wav というファイルを使います。このデータは3つの周波数が足し合わされたものでノイズも加えられているものです。サンプリング周波数は22050Hzで16ビットの4秒ほどのデータです。

f:id:terakei818:20160430042744p:plain

このファイルをwavread('sample.wav')で読み込みます(3行目)。この時sには波形データ、Fsにはサンプリング周波数、bitsにはビット数が格納されます。

次に音のデータの大きさを取得(5行目)←これがFFTのポイント数になる(みたい)

そしてFFTにはつきものの窓関数を作ります(8行目)。

窓関数というのは大雑把に説明すると、FFTなどを計算するときに無限の領域に対して計算することはできないためある区間からある区間までを区切ってサンプルをとり、その領域が繰り返し続く波と考えて計算をする。そのときの波がうまく周期的に分けられればいいが、もともと周波数を解析しようというのにうまく周期的に分けることはできないので両端がうまくつながらなくなってしまい(不連続になる)、フーリエ変換した結

果の周波数が歪んでしまうという困ったことが起きてしまう。

その問題を小さくするために区切った両端をうまくつないでくれる関数が窓関数であるのです。
詳しくは窓関数で調べてみてください。

ここでは解析によく使われる窓関数の一つのハミング窓というものを使っています。
窓関数を使うのは簡単で調べたい波形データ(ここではs)に窓関数を掛け合わせるだけです。

そして波形データをFFTをします(11行目)。

scilabでは便利な関数が用意されていてfft()という関数1つでFFTをしてくれます。めっちゃ簡単ですね。

今回は振幅特性を見るのでFFTした結果の絶対値をとります(13行目)。

これでグラフにプロットして終わりかと思いました、はじめは。

しかし世の中そんなに甘くはなかったです。このままグラフにプロットするとこんな感じになってしまいます。

f:id:terakei818:20160430032152p:plain

もうわけわかんないですね。サンプリング周波数は22050Hzなのに90000Hzまで出ていますし、3つの周波数の組み合わせのはずが6つもあります。なんか左右対称の雰囲気がありますね。←重要

お気づきかもしれませんがこのグラフ横のスケールが合っていません。なので今から横のスケールを周波数にします。

ここで重要となるのは周波数分解能というものです。

これはどれだけ細かい周波数まで見れるかというもので Fs/N で表されるものです。

Fsはサンプリング周波数、NはFFTのポイント数と言われFFTした後のサンプル数です。このNの値が大きいほど分解の性能は上がります。

グラフの1メモリが分解能にあたるわけですね。

またFFT(DFT)の性質として振幅特性ではN/2で左右対称になるというものがあります。先ほど失敗したグラフを見てもらえれば確かに左右対称になっていますね。

なので使うデータはN/2まででいいのです。

これをコードで書くと17行目のようになります。

後はグラフにプロットしてやればこんな感じになります。

f:id:terakei818:20160430041427p:plain

横軸のスケールが合っていい感じになりましたね。右端がナイキスト周波数までになっているのがばっちりです。

ナイキスト周波数とは周波数解析できる限界の周波数でサンプリング周波数の半分の値になります。

これにより含まれていた周波数は494Hz, 587Hz, 784Hzであると推測できました。(あくまで推測)

ここまで理解するのはなかなか大変でした。特にスケールを合わせるところはきちっとした説明が書いていないサイトが多くて苦労しました(;´д｀)トホホ

今回は長くなりましたのでいったんここで終わります。後編では実際にフィルタを作るところをやりたいと思います。