/*---------------------------------------------
*プログラム名 ： パワースペクトル
*ファイル名 ： ffttest.c
*バージョン ： up
*作成日 ： 1994/4
*-------------------------------------------- */
#include <stdio.h>
#include <process.h>
#include <conio.h>
#include <ctype.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323826
#define N 256 /* データの個数 */
int fft( int, double *, double * );
int main( void )
{
int j,i, n=256;
static double x[N], y[N],P[N]; /* data */
FILE *fp1;
for(j=0;j<256;j++){
x[j]=sin(2*PI*5*j*0.01)+sin(2*PI*20*j*0.01)+sin(2*PI*35*j*0.01);
}
for(j=0;j<256;j++)
y[j]=0;
fft( n, x, y );
for(i=1;i<=128;i++){
P[i]=0.00390625*(x[i]*x[i]+y[i]*y[i]);
printf("%f",y[i]);
putchar('\n');
}
fp1=fopen("fu-rie","w");
fprintf(fp1,"fu-rie \n");
for(i=1;i<=128;i++){
fprintf(fp1,"%3d%15f \n",i,P[i]);
fprintf(fp1,"\n");
}
fclose(fp1);
return( 0 );
}
/*-----------------------------------------------------------------------
* 関数名: int fft( int n, double x[], double y[] )
* 機 能: ＦＦＴ計算
* 入 力: n データ数
* : x データ配列（実部）
* : y データ配列（虚部）
* 出 力: 0
* 作成日: 1991/09/24
*---------------------------------------------------TextProgram--------*/
static void make_sintbl( int n, double sintbl[] )
{
int i, n2, n4, n8;
double c, s, dc, ds, t;
n2 = n / 2;
n4 = n / 4;
n8 = n / 8;
t = sin( PI / n );
dc = 2 * t * t;
ds = sqrt( dc * (2-dc) );
t = 2 * dc;
c = sintbl[n4] = 1.0;
s = sintbl[0] = .0;
for ( i = 1; i < n8; i++) {
c -= dc;
dc += t * c;
s += ds;
ds -= t * s;
sintbl[i] = s;
sintbl[n4-i] = c;
}
if ( n8 != 0 )
sintbl[n8] = sqrt(0.5);
for ( i = 0; i < n4; i++ )
sintbl[n2-i] = sintbl[i];
for ( i = 0; i < n2 + n4; i++ )
sintbl[i+n2] = -sintbl[i];
}
static void make_bitrev( int n, int bitrev[] )
{
int i, j, k, n2;
n2 = n / 2;
i = j = 0;
for ( ; ; ) {
bitrev[i] = j;
if ( ++i >= n )
break;
k = n2;
while ( k <= j ) {
j -= k;
k /= 2;
}
j += k;
}
}
int fft( int n, double x[], double y[] )
{
static int last_n = 0; /* 前回呼出時のn */
static int *bitrev = NULL; /* ビット反転表 */
static double *sintbl = NULL; /* 三角関数表 */
int i, j, k, ik, h, d, k2, n4, inverse;
double t, s, c, dx, dy;
/* 準備 */
inverse = 0;
n4 = n / 4;
if ( n != last_n || n == 0 ) {
last_n = n;
if ( sintbl != NULL )
free(sintbl);
if ( bitrev != NULL )
free(bitrev);
if ( n == 0 )
return 0; /* 記憶領域を解放した */
sintbl = malloc( (n+n4) * sizeof(double) );
bitrev = malloc( n * sizeof(int) );
if ( sintbl == NULL || bitrev == NULL ) {
fprintf( stderr, "記憶領域不足\n" ); return 1;
}
make_sintbl( n, sintbl );
make_bitrev( n, bitrev );
}
for ( i = 0; i < n; i++ ) { /* ビット反転 */
j = bitrev[i];
if ( i < j ) {
t = x[i]; x[i] = x[j]; x[j] = t;
t = y[i]; y[i] = y[j]; y[j] = t;
}
}
for ( k = 1; k < n; k = k2 ) { /* 変換 */
h = 0; k2 = k + k; d = n / k2;
for ( j = 0; j < k; j++ ) {
c = sintbl[h+n4];
s = sintbl[h];
for ( i = j; i < n; i += k2 ) {
ik = i + k;
dx = s * y[ik] + c * x[ik];
dy = c * y[ik] - s * x[ik];
x[ik] = x[i] - dx; x[i] += dx;
y[ik] = y[i] - dy; y[i] += dy;
}
h += d;
}
}
return( 0 );
}
の
fp1=fopen("fu-rie","w");
fprintf(fp1,"fu-rie \n");
for(i=1;i<=128;i++){
fprintf(fp1,"%3d%15f \n",i,P[i]);の
4行目のiの位置に周波数がくるらしいのですが、
周波数がどうだせばいいかわかりません。for分で128回繰り返すらしいのですが教えてください。
>周波数がどうだせばいいかわかりません
これってプログラム云々の問題ではないのではないですか？
周波数を出すために必要なものは時間と振動数ですよね。
単位時間の振動数が周波数なわけですから・・・
おそらく画像処理のプログラムかなにかを作ろうとしているのだと思いますが、おっしゃっている「周波数の出し方」というのが、私の認識では「平面画像」に関する事をさしているのであれば、「デジタル画像処理の専門書などを読んでください」が回答になってしまいますね・・・
ただソースを出して、「周波数がどうだせばいいかわかりません」ではなく、プログラムとしてどの部分でなにができないかを明確にしないと、誰も振り向いてくれないと思いますよ〜。
周波数をiの部分に入れるのですが、周波数は⊿ｆで、
プリントは⊿ｆ間隔で128回分 上限が128⊿ｆのグラフを書かないといけないらしいのです。
縦軸がパワースペクトル、横軸が周波数のグラフを書かないといけないのです。
⊿ｆ＝1/Ｎ⊿ｔ Ｎ＝256 ⊿ｔ＝0.01
です。
…意味解ってないでしょ。
256個のデータを処理してますが、これが何秒分のデータなのかが
わからんと周波数を算出できんってことが解ってないでしょ。
ちょっと厳しい言い方かもしれませんが・・・
まずはFFTの中身を分かっていますか？
FFTの中身がわかっているのであれば、「周波数を出す」ではなく、どのよう
な値からなにを導き出したいかを書き込んでもらわないと回答できないです。
もし、FFTの中身を理解していないのであれば、質問する場所が違います。
まずは画像処理や数学に関する書籍やHPなどで勉強してください。
その上で、アルゴリズムを考え、ソースを作成してどこがわからないのかを
明確にしてから、もう一度書き込んだ方が良いと思います。
そもそも、このソースは自分で作成したものですか？
自分で作成したものでないのであれば、まずこのソースとFFTのアルゴリズム
を理解した上で質問するべきですよね？
質問する本人が内容を理解していないものに対して質問をしても、回答する
側はなにを答えていいのかわからないですよ。