FFTの値がおかしくなる

Question

解決したいこと

バイナリファイルからデータを読み込んで配列に格納し、それをFFTにかけた後、スペクトル表示するというプログラムを書いていて、コンパイルして実行することは出来るのですが、FFTの結果がおかしくなってしまいます。
調べて色々試したのですが、解決しません。何が原因なのでしょうか。ご教授のほど、よろしくお願いいたします。
バイナリファイルの中身は以下のような数値です。
0.008271
-0.001495
0.000000
-0.007813
-0.014161
-0.026368
0.000000
0.037080
-0.010254
0.000488
-0.002960
・
・
・

発生している問題・エラーメッセージ

206099002025443869085040175710529487143113416799669242265250740530784206949264957106047804087609843367675182983333322436823044954917591898531079709544939520.000000 16174789132392824822629506959694874941668552674778128904007897595244676099903250209526602079076813033589364284034981126637158753898319687526134759555072.000000 110304749478410871324762407735844280167146792678333686310776179381347875861651481048045224304453279958148318539472420377016460201809005504620375874469888.000000 7393246484748170072092283098161425271834040984546707136614483038199951443564626579930399288747456525017394104736178698493855020572584283235181781877656649728.000000 
0.000000 
-0.067843
0.000000 
-0.008301 
-0.008301
-0.014649
0.000000 
-0.035127
-0.014649
-0.001007
-0.035127
-0.037080
0.000000 
-0.070742 
0.000000 
0.000977 
0.000000 
0.061983 
0.000000 
　　・
　　・
　　・

スペクトル表示

-90.275693
-86.033339
-86.722261
-87.587345
-91.649984
-85.913491
-86.367346
-87.104887
-89.295487
-88.738500
-89.742001
-87.371361
-88.316681
-88.074504
-87.900553
-92.438832
2960.794929
2962.997791
2962.997791
2960.794928
2960.794929
2962.997792
2962.997791
2960.794928
2960.794929
2962.997791
2962.997791
2960.794928
2960.794929
2962.997792
2962.997792
2960.794928

サンプリング点数が32の時の結果です。17個目の値からおかしくなってしまいます。

該当のソースコード

#define _CRT_SECURE_NO_WARNING
#include	<stdio.h>
#include	<stdlib.h>
#include     <string.h>
#include	<math.h>

int fft1();
int window_func();

int main(int argc,const char *argv[])
{
  int  frame; 　// 繰り返し回数
  int Nf = 1024; 　//サンプリング点数
  
  if(argc != 4)
  {
    printf("Usage: ./fft2 <No.1 input file name> <number of frames>\n");
    exit(1);
  }
 
  frame = atoi(argv[3]);
    
  int  i, j, k, n, iter, flag;
  
  char  data[100];
        
  double  ar1[Nf], ai1[Nf];
 
  double  tempr1[Nf], tempi1[Nf], tempr2[Nf], tempi2[Nf];
  
  double  dcnt[Nf];
  
  double  S1r1[Nf], S2r1[Nf], S3r1[Nf], S4r1[Nf];
  double  S1i1[Nf], S2i1[Nf], S3i1[Nf], S4i1[Nf];
  double  Myu1r1[Nf], Myu2r1[Nf], Myu3r1[Nf], Myu4r1[Nf];
  double  Myu1i1[Nf], Myu2i1[Nf], Myu3i1[Nf], Myu4i1[Nf];
  
  double  modi_S2[Nf], modi_av[Nf], modi_power[Nf];
  	
  FILE  *fp, *fp1, *fp2;
    
   /* check the number of input argument */  
 
  // initialize arrays
  
  for(i = 0; i < Nf; i++)
  {
    ar1[i] = 0.0;
    ai1[i] = 0.0;
    
    S1r1[i] = 0.0;
    S1i1[i] = 0.0;
    S2r1[i] = 0.0;
    S2i1[i] = 0.0;
    S3r1[i] = 0.0;
    S3i1[i] = 0.0;
    S4r1[i] = 0.0;
    S4i1[i] = 0.0;
    
    dcnt[i] = 0.0;    
    
    Myu1r1[i] = 0.0;
    Myu1i1[i] = 0.0;
    Myu2r1[i] = 0.0;
    Myu2i1[i] = 0.0;
    Myu3r1[i] = 0.0;
    Myu3i1[i] = 0.0;
    Myu4r1[i] = 0.0;
    Myu4i1[i] = 0.0;
     
    modi_S2[i] = 0.0;
    modi_av[i] = 0.0;
    
    modi_power[i] = 0.0;
  }  
  	
    /* read file for n_frame */

  fp1 = fopen(argv[1],"r");
  fp2 = fopen(argv[2],"r");
     
  iter = 0;
  flag = 0;
  
  for(n = 0; n < frame; n++)
  {
    for(i = 0; i < Nf; i++)
    {
      fread(&ar1[i], sizeof(double), 1, fp1);
    }

    for(i = 0; i < Nf; i++)
    {
      fread(&ai1[i], sizeof(double), 1, fp2);
    }
            
    window_func(ar1, Nf);　//窓関数
    window_func(ai1, Nf);　//窓関数
    fft1(ar1, ai1, Nf, iter, flag);　//FFT
    
    for(i = 0; i < Nf/2; i++)　//周波数分布の入れ替え
    {
      tempr1[i]  = ar1[i];
      ar1[i]     = ar1[Nf/2+i];
      ar1[Nf/2+i] = tempr1[i];
	     
      tempi1[i]  = ai1[i];
      ai1[i]     = ai1[Nf/2+i];
      ai1[Nf/2+i] = tempi1[i];
    }

    for(i = 0; i < Nf; i++)　//frameの数だけ繰り返し足す
    {
      S2r1[i]  +=  pow(ar1[i], 2.0);
      S2i1[i]  +=  pow(ai1[i], 2.0);

      dcnt[i]  +=  1.0;
    }
  }
  
  fclose(fp2);

  fclose(fp1);
  
  printf("\n");
  
  for(i = 0; i < Nf; i++)
  {
    Myu2r1[i] = S2r1[i] / dcnt[i];　//平均を取る
    Myu2i1[i] = S2i1[i] / dcnt[i];
  }

  for(i = 0; i < Nf; i++)
    modi_S2[i] = Myu2r1[i] + Myu2i1[i];
  
  for(i = 0; i < Nf; i++)　//スペクトル表示（対数で）
  {
    modi_av[i] = modi_S2[i] / 1000000;
    modi_power[i] = 10 * log10(modi_av[i]);
  }
  
  fp = fopen("freq.dat", "w"); //ファイルに出力
       
  for(i = 0; i < Nf; i++)
    fprintf(fp, "%lf\n", modi_power[i]);
    
  fclose(fp);
  
  return 0;
	
}

int fft1(ar,ai,n,iter,flag)　//FFT
  float ar[],ai[];
  int   n, iter, flag;
{
  int    i, j, it, xp, xp2, k, j1, j2, im1, jm1;
  double sign, w, wr, wi, dr1, dr2, di1, di2, tr, ti, arg;

  if(n < 2) return(999);

  if(iter <= 0)
  {
    iter = 0;
    i = n;
    
    while(1)
    {
      if((i /= 2) == 0) break;
      iter++;
    }
  }

  j = 1;
  
  for(i = 0; i < iter; i++) j *= 2;
  
  if(n != j) return(1);
  
  if(flag == 1) sign=1.0;
  else sign=-1.0;
  
  xp2 = n;

  for(it = 0; it < iter; it++)
  {
    xp = xp2;
    xp2 = xp / 2;
    w = 3.141592653589793 / xp2;
    
    for(k = 0; k < xp2; k++)
    {
      arg = k * w;
      wr = cos(arg);
      wi = sign * sin(arg);
      i = k - xp;
      
      for(j = xp; j <= n; j += xp)
      {
        j1  = j + i;
        j2  = j1 + xp2;
        dr1 = ar[j1];
        dr2 = ar[j2];
        di1 = ai[j1];
        di2 = ai[j2];
        tr  = dr1 - dr2;
        ti  = di1 - di2;
        
        ar[j1] = dr1 + dr2;
        ai[j1] = di1 + di2;
        ar[j2] = tr * wr - ti * wi;
        ai[j2] = ti * wr + tr * wi;	      
      }
    }
  }

  j1 = n / 2;
  j2 = n - 1;
  j = 1;

  for(i = 1; i <= j2; i++)
  {
    if(i < j)
    {
      im1 = i - 1;
      jm1 = j - 1;
      tr  = ar[jm1];
      ti  = ai[jm1];
      
      ar[jm1] = ar[im1];
      ai[jm1] = ai[im1];
      ar[im1] = tr;
      ai[im1] = ti;
    }
  
    k = j1;
  
    while(k < j)
    {
      j -= k;
      k /= 2;
    }
	  
    j += k;
  }

  if(flag == 0) return(0);
  
  w = n;

  for(i = 0; i < n; i++)
  {
    ar[i] = ar[i] / w;
    ai[i] = ai[i] / w;
  }
  
  return(0);
}


int window_func(frame,n)　//窓関数
  float frame[];
  int	n;
{
  int	 i;
  float  winv[n];
  double pi = 3.141592653589793;

  for(i = 0; i < n; i++)
  {
    /* Blackman
    winv[i] = 0.42 - 0.5 * cos(2.0*pi*i/(n-1)) + 0.08 * cos(4.0*pi*i/(n-1)); */

    /* Blackman-Harris */
    winv[i] = 0.35875 - 0.48829 * cos(2.0*pi*i/(n-1)) + 0.14128 * cos(4.0*pi*i/(n-1)) - 0.01168 * cos(6.0*pi*i/(n-1));

    /* no window*/
    /*winv[i]=1.0;*/ 

    frame[i] = frame[i] * winv[i];
  }

  return(0);
}

Comments

[thkana]

> 結果がおかしくなってしまいます。
「おかしい」という判断はどうやってしましたか。

> 調べて色々試したのですが
頑張っているアピールをするためではなく、解決に向けて有効な情報とするために
何を試した->どういう結果を得た->結果からどういうことがわかった(ここはおかしくない/ここは少なくともおかしい)、
などということを記載してください。色々試したのなら、少なくともやったことと得た結果は記載できるはずです。

[jimbe]

>調べて色々試したのですが、解決しません。何が原因なのでしょうか。

teratail は依頼をするところではありません。
どこまで正常に計算されていることを確認されましたか。
17個目からおかしいのであれば、17個目の時の各変数の値が想定通りなのかは確認されたのでしょうか。

…とりあえずコードを提示するならコンパイルエラーは出ない状態の方が良いと思います。(データファイルが無ければ実行しても再現出来ませんが…)

[yu_89]

thkana様、jimbe様、ご忠告いただき、ありがとうございます。
どこでエラーとなっているのか分からず、丸投げな質問になってしまい、申し訳ございませんでした。
次回からはもっと詳細に記述するように気を付けます。
確かに、コードがあってもデータファイルがファイルがなければ再現できませんよね。それも提示するのを忘れていました。コードも貼り付ける時にコンパイルエラーが出ないことを確認してから提示します。
色々書くことが足りず、本当に申し訳ございませんした。

Answer 1

fftと窓関数の変数の型をfloatからdoubleに書き換えたら直りました。
お騒がせして申し訳ございませんでした。

Comments

[thkana]

「今の」Cの文法で書きましょうよ...そうしていればそもそも防げた(検出できた)間違いです。
int fft1(ar,ai,n,iter,flag)　//FFT
float ar[],ai[];
int n, iter, flag;
{
これ、1989年版規格以前の記法です。

[yu_89]

ご回答いただき、ありがとうございます。
そうなんですね！それは知りませんでした。本に載ってたコードをそのまま使ってました。
そうですね、今の見慣れた文法で書いていれば、見つけられる間違いだと思います。すぐに書き直します。
教えていただき、ありがとうございます！