vDSPを使う(高速フーリエ変換編 その3)


前の記事で、高速フーリエ変換を行うことができました。とはいえ、フーリエ変換された結果は複素数の配列になっているので、周波数特性や位相特性を見るためにはもう一段階処理をする必要があります。

まずはデータを観察してみる

512個も要素があるので少し見にくいですが、とりあえずフーリエ変換した結果がどうなっているのか見てみましょう。(当然、これはsplitComplexの中身をfreeする前に行います)

// 結果を表示してみる
for (int i=0; i<512; i++) {
    NSLog(@"[%d] %.2f, %.2f",i, splitComplex.realp[i],
        splitComplex.imagp[i]);
}

すると、大抵の所は0に近い値が並んでいて、30, 60, 452, 482の所だけ数値が飛び出ていることがわかります。おや、この30とか60という数字、ちょっと前にどこかで見覚えがありませんか?

ログを見てるだとわかりにくいですが、482 = 512 – 30、452 = 512 – 60なのでこの配列は真ん中を中心に、実数(real)が線対称、虚数(imag)が点対称になっています。

面倒くさいので結論から言うと、半分より向こう側(256〜511)は使いません

で、この配列は?

音声信号をフーリエ変換した結果取得できるのは、周波数特性と位相特性という二つのものです。「なるほど、ではrealかimagのどっちかが周波数特性で、もう一つが位相特性になっているのだな」と、そうなっていればよいのですが、残念ながらもう一段階変換を行わないと周波数特性と位相特性は得られません。

複素数の配列から複素平面を思い浮かべたできた人なら予想できたかもしれませんが、振幅から周波数特性、角度から位相特性を取得することができます。

周波数特性

とりあえず今回はこれから周波数特性を取得してみます。位相特性ってなんか地味だし。

上にも書いた通り、複素数の振幅=その周波数がどれくらい含まれているかを表すので、その複素数の振幅を求めてやります。(ちなみに、ここの処理もvDSP_vdistという関数を使うと高速にできるので、挑戦してみてください。ここではわかりやすさ重視で普通に書きます。)

// 結果を表示してみる
for (int i=0; i<256; i++) {// 半分よりあっち側は無意味なので無視するよ!
	float real = splitComplex.realp[i];
	float imag = splitComplex.imagp[i];
	float distance = sqrt(real*real + imag*imag);
	NSLog(@"[%d] %.2f",i, distance);
}

すると、30と60のところにそれぞれ256, 128という数字がでています。波形を生成する部分で、30回振動するサイン波と半分の大きさでその倍(60回)振動するサイン波が加算されているので、その特性がここにあらわれています。

試しに、周波数を他の数字にしてみて、解析結果に反映されるか確かめてみましょう。

おや?

周波数が整数の時は比較的きれいに(本来含まれていない周波数は0になっている)解析できますが、中途半端な少数の時はあまりキレイに分離されません。

離散フーリエ変換ではそのアルゴリズムの特性上、どうしても信号の長さをキレイに割り切れる周波数以外の成分はきれいに解析できません。しかし、それを改善するための方法があるのですが、長くなったのでまた次の記事で。

参考リンク


“vDSPを使う(高速フーリエ変換編 その3)” への 3 件のフィードバック

  1. [0] [1] [2]の数字の部分の単位はHzでいいんですよね?

    計算する周波数を、もう少し細かく[1.1] [1.2] [1.3] [1.4]・・・と、細かく結果を計算したいのですが、[0] [1] [2] と整数の時の結果しか計算できないんでしょうか?

  2. 単位についてはサンプリングレートによって変動します。ここで書いているサンプルでは、サンプリングレート=512の場合はHzになります。

    基本的にフーリエ解析は倍音成分を分析するものなので、そのような小数倍で成分を出す事はできないです。
    サンプルの長さを長めに取ることで、解析できる倍音の数が増えるのでそのような方法はどうでしょう。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です