C言語1
1#include <stdio.h>2#include <stdlib.h>3#include <math.h>4 5int fft1();6int window_func();7double* Kurtosis(double a[], double b[], double c[], double d[], int e);8 9int main(int argc,const char *argv[])10{11 int N = 128;12 int i, j, k, n, n_frame, iter, flag;13 14 char datar[256], datai[256];15 16 float ar1[N], ai1[N];17 18 double tempr1[N], tempi1[N];19 20 double pwr1, pwi1;21 22 double dcnt[N];23 24 double S1r1[N], S2r1[N], S3r1[N], S4r1[N];25 double S1i1[N], S2i1[N], S3i1[N], S4i1[N];26 double Myu1r1[N], Myu2r1[N], Myu3r1[N], Myu4r1[N];27 double Myu1i1[N], Myu2i1[N], Myu3i1[N], Myu4i1[N];28 double M1r[N], M2r[N], M3r[N], M4r[N];29 double M1i[N], M2i[N], M3i[N], M4i[N]; 30 31 double* K = NULL;32 33 FILE *fp, *fp1, *fp2, *fp3;34 35 /* check the number of input argument */36 37 if(argc != 4)38 {39 printf("Usage: ./fft2 <No.1 input file name> <No.2 input file name> <number of frames>\n");40 exit(1);41 }42 43 n_frame=atoi(argv[3]);44 45 // initialize arrays46 47 for(i = 0; i < N; i++)48 {49 ar1[i] = 0.0;50 ai1[i] = 0.0;51 52 M1r[i] = 0.0;53 M1i[i] = 0.0;54 M2r[i] = 0.0;55 M2i[i] = 0.0;56 M3r[i] = 0.0;57 M3i[i] = 0.0;58 M4r[i] = 0.0;59 M4i[i] = 0.0;60 61 S1r1[i] = 0.0;62 S1i1[i] = 0.0;63 S2r1[i] = 0.0;64 S2i1[i] = 0.0;65 S3r1[i] = 0.0;66 S3i1[i] = 0.0;67 S4r1[i] = 0.0;68 S4i1[i] = 0.0;69 70 dcnt[i] = 0.0; 71 72 tempr1[i] = 0.0;73 tempi1[i] = 0.0;74 75 Myu1r1[i] = 0.0;76 Myu1i1[i] = 0.0;77 Myu2r1[i] = 0.0;78 Myu2i1[i] = 0.0;79 Myu3r1[i] = 0.0;80 Myu3i1[i] = 0.0;81 Myu4r1[i] = 0.0;82 Myu4i1[i] = 0.0;83 }84 85 pwr1 = 0.0;86 pwi1 = 0.0;87 88 /* read file for n_frame */89 90 fp1 = fopen(argv[1],"r");91 fp2 = fopen(argv[2],"r");92 93 iter = 0;94 flag = 0;95 96 for(n = 0; n < n_frame; n++)97 {98 printf("\rProcessing %10d th frame",n+1);99 100 101 for(i = 0; i < N; i++)102 {103 fgets(datar, sizeof(datar), fp1);104 ar1[i] = atof(datar);105 }106 107 for(i = 0; i < N; i++)108 {109 fgets(datai, sizeof(datai), fp2);110 ar1[i] = atof(datai);111 }112 113 printf("\n");114 115 window_func(ar1, N);116 window_func(ai1, N);117 fft1(ar1, ai1, N, iter, flag);118 119 for(i = 0; i < N/2; i++)120 {121 tempr1[i] = ar1[i];122 ar1[i] = ar1[N/2+i];123 ar1[N/2+i] = tempr1[i];124 125 tempi1[i] = ai1[i];126 ai1[i] = ai1[N/2+i];127 ai1[N/2+i] = tempi1[i];128 }129 130 for(i = 0; i < N; i++)131 {132 pwr1 = ar1[i];133 pwi1 = ai1[i];134 135 S1r1[i] += (double)pwr1;136 S1i1[i] += (double)pwi1;137 138 S2r1[i] += pow((double)pwr1, 2.0);139 S2i1[i] += pow((double)pwi1, 2.0);140 141 S3r1[i] += pow((double)pwr1, 3.0);142 S3i1[i] += pow((double)pwi1, 3.0);143 144 S4r1[i] += pow((double)pwr1, 4.0);145 S4i1[i] += pow((double)pwi1, 4.0);146 147 dcnt[i] += 1.0;148 }149 }150 151 printf("\n");152 fclose(fp1);153 fclose(fp2);154 155 for(i = 0; i < N; i++)156 {157 Myu1r1[i] = S1r1[i] / dcnt[i];158 Myu1i1[i] = S1i1[i] / dcnt[i];159 160 Myu2r1[i] = S2r1[i] / dcnt[i];161 Myu2i1[i] = S2i1[i] / dcnt[i];162 163 Myu3r1[i] = S3r1[i] / dcnt[i];164 Myu3i1[i] = S3i1[i] / dcnt[i];165 166 Myu4r1[i] = S4r1[i] / dcnt[i];167 Myu4i1[i] = S4i1[i] / dcnt[i];168 }169 170 K = Kurtosis(Myu1r1, Myu2r1, Myu3r1, Myu4r1, N); //エラーがでる箇所171 for(i = 0; i < N; i++)172 printf("%f\n", K[i]);173 174 free(K);175 176 return 0;177 178}179 180//FFTの計算181int fft1(ar,ai,n,iter,flag)182 float ar[],ai[];183 int n, iter, flag;184{185 int i, j, it, xp, xp2, k, j1, j2, im1, jm1;186 double sign, w, wr, wi, dr1, dr2, di1, di2, tr, ti, arg;187 188 if(n < 2) return(999);189 190 if(iter <= 0)191 {192 iter = 0;193 i = n;194 195 while(1)196 {197 if((i /= 2) == 0) break;198 iter++;199 }200 }201 202 j = 1;203 204 for(i = 0; i < iter; i++) j *= 2;205 206 if(n != j) return(1);207 208 if(flag == 1) sign=1.0;209 else sign=-1.0;210 211 xp2 = n;212 213 for(it = 0; it < iter; it++)214 {215 xp = xp2;216 xp2 = xp / 2;217 w = 3.141592653589793 / xp2;218 219 for(k = 0; k < xp2; k++)220 {221 arg = k * w;222 wr = cos(arg);223 wi = sign * sin(arg);224 i = k - xp;225 226 for(j = xp; j <= n; j += xp)227 {228 j1 = j + i;229 j2 = j1 + xp2;230 dr1 = ar[j1];231 dr2 = ar[j2];232 di1 = ai[j1];233 di2 = ai[j2];234 tr = dr1 - dr2;235 ti = di1 - di2;236 237 ar[j1] = dr1 + dr2;238 ai[j1] = di1 + di2;239 ar[j2] = tr * wr - ti * wi;240 ai[j2] = ti * wr + tr * wi; 241 }242 }243 }244 245 j1 = n / 2;246 j2 = n - 1;247 j = 1;248 249 for(i = 1; i <= j2; i++)250 {251 if(i < j)252 {253 im1 = i - 1;254 jm1 = j - 1;255 tr = ar[jm1];256 ti = ai[jm1];257 258 ar[jm1] = ar[im1];259 ai[jm1] = ai[im1];260 ar[im1] = tr;261 ai[im1] = ti;262 }263 264 k = j1;265 266 while(k < j)267 {268 j -= k;269 k /= 2;270 }271 272 j += k;273 }274 275 if(flag == 0) return(0);276 277 w = n;278 279 for(i = 0; i < n; i++)280 {281 ar[i] = ar[i] / w;282 ai[i] = ai[i] / w;283 }284 285 return(0);286}287 288//窓関数の処理を行う関数289int window_func(frame,n)290 float frame[];291 int n;292{293 int i;294 float winv[n];295 double pi = 3.141592653589793;296 297 for(i = 0; i < n; i++)298 {299 /* Blackman 300 winv[i] = 0.42 - 0.5 * cos(2.0*pi*i/(n-1)) + 0.08 * cos(4.0*pi*i/(n-1)); */301 302 /* Blackman-Harris */303 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));304 305 /* no window*/306 /*winv[i]=1.0;*/ 307 308 frame[i] = frame[i] * winv[i];309 }310 311 return(0);312}313 314//問題としている自作関数315double* Kurtosis(double Myu1[], double Myu2[], double Myu3[], double Myu4[], int N)316{317 int i;318 double M1[N], M2[N], M3[N], M4[N];319 double top[N], bottom[N];320 double* K = (double*)malloc(N);321 322 for(i = 0; i < N; i++)323 {324 M1[i] = 0.0;325 M2[i] = 0.0;326 M3[i] = 0.0;327 M4[i] = 0.0;328 329 top[i] = 0.0;330 bottom[i] = 0.0;331 }332 333 for(i = 0; i < N; i++)334 {335 M1[i] = 4.0 * Myu3[i] * Myu1[i];336 337 M2[i] = 6.0 * Myu2[i] * Myu1[i] * Myu1[i];338 339 M3[i] = 3.0 * Myu1[i] * Myu1[i] * Myu1[i] * Myu1[i];340 341 M4[i] = Myu2[i] - (Myu1[i] * Myu1[i]);342 }343 344 for(i = 0; i < N; i++)345 {346 top[i] = Myu4[i] - M1[i] + M2[i] - M3[i]; 347 bottom[i] = M4[i];348 }349 350 for(i = 0; i < N; i++)351 {352 K[i] = top[i] / bottom[i]; 353 }354}
0 コメント