S(k)=10*log{1/N||シグマ(上N-1下n=0)s(n)h(n)exp(-j2pink/N)||^2}
のプログラムを教えてください。これはパワースペクトルの式です。
N;512 h(n)；ハニング窓関数で(1/2{1-cos(2pin/N)})です。s(n)；オーディオ信号です。VisualC++でお願いします。今現在わかっているのは、シグマの中です。ここでa(j)というのは自分の式でいうh(n)s(n)をかけたものです。窓関数のプログラムもわかります。残りのものとどうひとつにまとめるかがわかりません。
#include <math.h>
void rfft(int n, double a[])
{
int m, mh, mq, i, j, k, jr, ji, kr, ki;
double theta, wr, wi, xr, xi;
i = 0;
for (j = 1; j < n - 1; j++) {
for (k = n >> 1; k > (i ^= k); k >>= 1);
if (j < i) {
xr = a[j];
a[j] = a[i];
a[i] = xr;
}
}
theta = -8 * atan(1.0); /* -2*pi */
for (mh = 1; (m = mh << 1) <= n; mh = m) {
mq = mh >> 1;
theta *= 0.5;
/* ---- real to real butterflies (W == 1) ---- */
for (jr = 0; jr < n; jr += m) {
kr = jr + mh;
xr = a[kr];
a[kr] = a[jr] - xr;
a[jr] += xr;
}
for (i = 1; i < mq; i++) {
wr = cos(theta * i);
wi = sin(theta * i);
for (j = 0; j < n; j += m) {
jr = j + i;
ji = j + mh - i;
kr = j + mh + i;
ki = j + m - i;
xr = wr * a[kr] + wi * a[ki];
xi = wr * a[ki] - wi * a[kr];
a[kr] = -a[ji] + xi;
a[ki] = a[ji] + xi;
a[ji] = a[jr] - xr;
a[jr] = a[jr] + xr;
}
}
タイトル「高速フーリエ変換からパワースペクトルを求める」を
見てください。
で、納得できたら解決をチェックしてね。
でも、コメントが足りませんね。
atan(1.0); が (1/4)πであることをサッとわかる人は少ないかもしれません。
# 私もこっち系です。3.14.... 書くの面倒ですよね。
# math.h にも PI の定義がなかったし。(VC++)
また、このようなアルゴリズム系は
すでにコーディングして公開されている方がいらっしゃるので、
そういったソースコードを読むことからはじめた方が
解決が早いのではないでしょうか。
普通 sin,cos はテーブル化しますよね。