プログラムを教えてください

質問

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)；オーディオ信号です。お願いします。

回答No.7

noname#5537

s(n) も自分で発生させるんですか。
オーディオ信号ということだったので，てっきりファイルなどから読み込むのかと思ってました。

> えっと・・ｓ（ｎ）＝sin(n)とかでも平気ですかね？

大丈夫でしょう。

補足

やっぱり出来ません。勉強不足かもしれませんが・・具体的にプログラムを教えてもらうことは無理ですか？その次からはもっと自分自身で勉強するんで・・２０日までにそのプログラムを作らないと卒業できなくなってしまうんです。その次の課題のプログラムは自分自身でしっかり勉強するんで・・お願いします。助けてください。って無理ですか？

回答No.6

noname#5537

> S(n)はオーディオ信号なのですが・・S（ｎ）はどうおけばよろしいでしょうか？

質問の意図が良く分かりません。
「どうおけば」とはどういう意味でしょうか？

補足

えっと・・ｓ（ｎ）＝sin(n)とかでも平気ですかね？

回答No.5

noname#5537

http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/ftmn2_12.html

これですか？
丸ごとコピーして掲載するのはまずいと思いますよ。

F[k]=Σ_j=0^n-1 a[j]*exp(-2*pi*i*j*k/n),0<=k<n …(1)
S(k)=10*log{1/N||シグマ(上N-1下n=0)s(n)h(n)exp(-j2pink/N)||^2} …(2)

(1)式の j は，(2)式の n に相当します。
(1)式の a(j) は (2)式の s(n)h(n) です。
(1)式の i は虚数単位です。

補足

ありがとうございます。S(n)はオーディオ信号なのですが・・S（ｎ）はどうおけばよろしいでしょうか？

回答No.4

noname#5537

はじめに訂正です。

> 最後に 10*log を掛けて dB (デシベル)単位に変換してます。

「log (常用対数) をとって 10 を掛ける」の間違いでした。

> 複素数の計算は実部と虚部にわけるのですか？？

基本的にはそういうことになります。

# 最新の C の規格では組み込み型として複素数が使えるようですが。
# また，C++ なら std::complex が使えます。(こちらは組み込み型ではないです)

補足

F[k]=Σ_j=0^n-1 a[j]*exp(-2*pi*i*j*k/n),0<=k<n を計算する.
出力 a[0...n/2], a[n/2+1...n-1] はそれぞれ F[0...n/2] の実部,
F[n/2+1...n-1] の虚部に相当する．
*/
#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;

/* ---- scrambler ---- */
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;
}
/* ---- complex to complex butterflies (W != 1) ---- */
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;
}
}
/* ---- real to complex butterflies are trivial ---- */
}

このプログラムを使おうと思うのですが・・１行目のシグマのあとのjと式の中のjはちがいますか？たびたび申しわありません。a(j)がh(n)s(n)にあたるんですよね？すいません。お願いします。

回答No.3

noname#5537

genkitake さん，この式の意味はわかっていますか？

まず離散フーリエ変換(DFT)を計算して，
その絶対値の二乗を計算して，
最後に 10*log を掛けて dB (デシベル)単位に変換してます。

DFT の計算は通常 FFT というアルゴリズムを使います，
勉強のためなら自分で組んでもいいですが，
FFT を計算してくれるライブラリが沢山ありますから，
そういうのを使うといいでしょう。
(参考 URL にもあります。)

窓関数は，最初に入力信号 s(n) に掛け合わせておけばいいだけです。

また何か分からないことがあったら質問してください。
ではがんばってください。

# ちなみに，j は虚数単位。複素数の計算になります。

参考URL：

http://momonga.t.u-tokyo.ac.jp/~ooura/fftman/

補足

ありがとうございます。プログラム初心者なもので・・・複素数の計算は実部と虚部にわけるのですか？？

回答No.2

επιστημη（@episteme）

ちょっとマチガイ。

[1] Σの内部 f(n): s(n)h(n)exp(-j2pink/N) を求める関数 double f(int n) を書きなさい。
[2] n = 0 .. N-1 に対し f(n) を求め、その結果を積算しなさい。
[3] [2]の結果をNで割り、二乗し、10を底とする対数を取って10倍すれば答が求まります。

お礼

ありがとうございます。がんばってみます。

回答No.1

επιστημη（@episteme）

[1] Σの内部 f(n): s(n)h(n)exp(-j2pink/N)||^2 を求める関数 double f(int n) を書きなさい。
[2] n = 0 .. N-1 に対し f(n) を求め、その結果を積算しなさい。
[3] [2]の結果をNで割り、10を底とする対数を取って10倍すれば答が求まります。

ここはあなたのかわりに誰かがプログラムを書いてくれるところではありません。