質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.69%
C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

Q&A

解決済

1回答

557閲覧

FFTの値がおかしくなる

yu_89

総合スコア34

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

0グッド

0クリップ

投稿2023/01/11 03:23

解決したいこと

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


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

FFTをかけた後

1206099002025443869085040175710529487143113416799669242265250740530784206949264957106047804087609843367675182983333322436823044954917591898531079709544939520.000000 16174789132392824822629506959694874941668552674778128904007897595244676099903250209526602079076813033589364284034981126637158753898319687526134759555072.000000 110304749478410871324762407735844280167146792678333686310776179381347875861651481048045224304453279958148318539472420377016460201809005504620375874469888.000000 7393246484748170072092283098161425271834040984546707136614483038199951443564626579930399288747456525017394104736178698493855020572584283235181781877656649728.000000 20.000000 3-0.067843 40.000000 5-0.008301 6-0.008301 7-0.014649 80.000000 9-0.035127 10-0.014649 11-0.001007 12-0.035127 13-0.037080 140.000000 15-0.070742 160.000000 170.000977 180.000000 190.061983 200.000000 21  ・ 22  ・ 23  ・

スペクトル表示

-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個目の値からおかしくなってしまいます。

該当のソースコード

fft.c

1#define _CRT_SECURE_NO_WARNING 2#include <stdio.h> 3#include <stdlib.h> 4#include <string.h> 5#include <math.h> 6 7int fft1(); 8int window_func(); 9 10int main(int argc,const char *argv[]) 11{ 12 int frame;  // 繰り返し回数 13 int Nf = 1024;  //サンプリング点数 14 15 if(argc != 4) 16 { 17 printf("Usage: ./fft2 <No.1 input file name> <number of frames>\n"); 18 exit(1); 19 } 20 21 frame = atoi(argv[3]); 22 23 int i, j, k, n, iter, flag; 24 25 char data[100]; 26 27 double ar1[Nf], ai1[Nf]; 28 29 double tempr1[Nf], tempi1[Nf], tempr2[Nf], tempi2[Nf]; 30 31 double dcnt[Nf]; 32 33 double S1r1[Nf], S2r1[Nf], S3r1[Nf], S4r1[Nf]; 34 double S1i1[Nf], S2i1[Nf], S3i1[Nf], S4i1[Nf]; 35 double Myu1r1[Nf], Myu2r1[Nf], Myu3r1[Nf], Myu4r1[Nf]; 36 double Myu1i1[Nf], Myu2i1[Nf], Myu3i1[Nf], Myu4i1[Nf]; 37 38 double modi_S2[Nf], modi_av[Nf], modi_power[Nf]; 39 40 FILE *fp, *fp1, *fp2; 41 42 /* check the number of input argument */ 43 44 // initialize arrays 45 46 for(i = 0; i < Nf; i++) 47 { 48 ar1[i] = 0.0; 49 ai1[i] = 0.0; 50 51 S1r1[i] = 0.0; 52 S1i1[i] = 0.0; 53 S2r1[i] = 0.0; 54 S2i1[i] = 0.0; 55 S3r1[i] = 0.0; 56 S3i1[i] = 0.0; 57 S4r1[i] = 0.0; 58 S4i1[i] = 0.0; 59 60 dcnt[i] = 0.0; 61 62 Myu1r1[i] = 0.0; 63 Myu1i1[i] = 0.0; 64 Myu2r1[i] = 0.0; 65 Myu2i1[i] = 0.0; 66 Myu3r1[i] = 0.0; 67 Myu3i1[i] = 0.0; 68 Myu4r1[i] = 0.0; 69 Myu4i1[i] = 0.0; 70 71 modi_S2[i] = 0.0; 72 modi_av[i] = 0.0; 73 74 modi_power[i] = 0.0; 75 } 76 77 /* read file for n_frame */ 78 79 fp1 = fopen(argv[1],"r"); 80 fp2 = fopen(argv[2],"r"); 81 82 iter = 0; 83 flag = 0; 84 85 for(n = 0; n < frame; n++) 86 { 87 for(i = 0; i < Nf; i++) 88 { 89 fread(&ar1[i], sizeof(double), 1, fp1); 90 } 91 92 for(i = 0; i < Nf; i++) 93 { 94 fread(&ai1[i], sizeof(double), 1, fp2); 95 } 96 97 window_func(ar1, Nf); //窓関数 98 window_func(ai1, Nf); //窓関数 99 fft1(ar1, ai1, Nf, iter, flag); //FFT 100 101 for(i = 0; i < Nf/2; i++) //周波数分布の入れ替え 102 { 103 tempr1[i] = ar1[i]; 104 ar1[i] = ar1[Nf/2+i]; 105 ar1[Nf/2+i] = tempr1[i]; 106 107 tempi1[i] = ai1[i]; 108 ai1[i] = ai1[Nf/2+i]; 109 ai1[Nf/2+i] = tempi1[i]; 110 } 111 112 for(i = 0; i < Nf; i++) //frameの数だけ繰り返し足す 113 { 114 S2r1[i] += pow(ar1[i], 2.0); 115 S2i1[i] += pow(ai1[i], 2.0); 116 117 dcnt[i] += 1.0; 118 } 119 } 120 121 fclose(fp2); 122 123 fclose(fp1); 124 125 printf("\n"); 126 127 for(i = 0; i < Nf; i++) 128 { 129 Myu2r1[i] = S2r1[i] / dcnt[i]; //平均を取る 130 Myu2i1[i] = S2i1[i] / dcnt[i]; 131 } 132 133 for(i = 0; i < Nf; i++) 134 modi_S2[i] = Myu2r1[i] + Myu2i1[i]; 135 136 for(i = 0; i < Nf; i++) //スペクトル表示(対数で) 137 { 138 modi_av[i] = modi_S2[i] / 1000000; 139 modi_power[i] = 10 * log10(modi_av[i]); 140 } 141 142 fp = fopen("freq.dat", "w"); //ファイルに出力 143 144 for(i = 0; i < Nf; i++) 145 fprintf(fp, "%lf\n", modi_power[i]); 146 147 fclose(fp); 148 149 return 0; 150 151} 152 153int fft1(ar,ai,n,iter,flag) //FFT 154 float ar[],ai[]; 155 int n, iter, flag; 156{ 157 int i, j, it, xp, xp2, k, j1, j2, im1, jm1; 158 double sign, w, wr, wi, dr1, dr2, di1, di2, tr, ti, arg; 159 160 if(n < 2) return(999); 161 162 if(iter <= 0) 163 { 164 iter = 0; 165 i = n; 166 167 while(1) 168 { 169 if((i /= 2) == 0) break; 170 iter++; 171 } 172 } 173 174 j = 1; 175 176 for(i = 0; i < iter; i++) j *= 2; 177 178 if(n != j) return(1); 179 180 if(flag == 1) sign=1.0; 181 else sign=-1.0; 182 183 xp2 = n; 184 185 for(it = 0; it < iter; it++) 186 { 187 xp = xp2; 188 xp2 = xp / 2; 189 w = 3.141592653589793 / xp2; 190 191 for(k = 0; k < xp2; k++) 192 { 193 arg = k * w; 194 wr = cos(arg); 195 wi = sign * sin(arg); 196 i = k - xp; 197 198 for(j = xp; j <= n; j += xp) 199 { 200 j1 = j + i; 201 j2 = j1 + xp2; 202 dr1 = ar[j1]; 203 dr2 = ar[j2]; 204 di1 = ai[j1]; 205 di2 = ai[j2]; 206 tr = dr1 - dr2; 207 ti = di1 - di2; 208 209 ar[j1] = dr1 + dr2; 210 ai[j1] = di1 + di2; 211 ar[j2] = tr * wr - ti * wi; 212 ai[j2] = ti * wr + tr * wi; 213 } 214 } 215 } 216 217 j1 = n / 2; 218 j2 = n - 1; 219 j = 1; 220 221 for(i = 1; i <= j2; i++) 222 { 223 if(i < j) 224 { 225 im1 = i - 1; 226 jm1 = j - 1; 227 tr = ar[jm1]; 228 ti = ai[jm1]; 229 230 ar[jm1] = ar[im1]; 231 ai[jm1] = ai[im1]; 232 ar[im1] = tr; 233 ai[im1] = ti; 234 } 235 236 k = j1; 237 238 while(k < j) 239 { 240 j -= k; 241 k /= 2; 242 } 243 244 j += k; 245 } 246 247 if(flag == 0) return(0); 248 249 w = n; 250 251 for(i = 0; i < n; i++) 252 { 253 ar[i] = ar[i] / w; 254 ai[i] = ai[i] / w; 255 } 256 257 return(0); 258} 259 260 261int window_func(frame,n) //窓関数 262 float frame[]; 263 int n; 264{ 265 int i; 266 float winv[n]; 267 double pi = 3.141592653589793; 268 269 for(i = 0; i < n; i++) 270 { 271 /* Blackman 272 winv[i] = 0.42 - 0.5 * cos(2.0*pi*i/(n-1)) + 0.08 * cos(4.0*pi*i/(n-1)); */ 273 274 /* Blackman-Harris */ 275 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)); 276 277 /* no window*/ 278 /*winv[i]=1.0;*/ 279 280 frame[i] = frame[i] * winv[i]; 281 } 282 283 return(0); 284}

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

下記のような質問は推奨されていません。

  • 質問になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

適切な質問に修正を依頼しましょう。

また依頼した内容が修正された場合は、修正依頼を取り消すようにしましょう。

thkana

2023/01/11 04:37

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

2023/01/11 07:12 編集

>調べて色々試したのですが、解決しません。何が原因なのでしょうか。 teratail は依頼をするところではありません。 どこまで正常に計算されていることを確認されましたか。 17個目からおかしいのであれば、17個目の時の各変数の値が想定通りなのかは確認されたのでしょうか。 …とりあえずコードを提示するならコンパイルエラーは出ない状態の方が良いと思います。(データファイルが無ければ実行しても再現出来ませんが…)
yu_89

2023/01/11 07:36

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

回答1

0

自己解決

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

投稿2023/01/11 07:38

yu_89

総合スコア34

下記のような回答は推奨されていません。

  • 質問の回答になっていない投稿
  • スパムや攻撃的な表現を用いた投稿

このような回答には修正を依頼しましょう。

また依頼した内容が修正された場合は、修正依頼を取り消すようにしましょう。

回答へのコメント

thkana

2023/01/11 11:44

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

2023/01/12 03:30

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.69%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

同じタグがついた質問を見る

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。