test.c
#include <stdio.h> #include <stdlib.h> #include <math.h> int main(int argc,const char *argv[]) { int N=1024; int i, j, n, n_frame, iter, flag; float ar1[N], ai1[N], ar2[N], ai2[N]; double tempr1[N], tempi1[N], tempr2[N], tempi2[N]; double pwr1, pwi1, pwr2, pwi2; double dcnt[N]; double S1r1[N], S2r1[N], S3r1[N], S4r1[N]; double S1i1[N], S2i1[N], S3i1[N], S4i1[N]; double Myu1r1[N], Myu2r1[N], Myu3r1[N], Myu4r1[N]; double Myu1i1[N], Myu2i1[N], Myu3i1[N], Myu4i1[N]; double S1r2[N], S2r2[N], S3r2[N], S4r2[N]; double S1i2[N], S2i2[N], S3i2[N], S4i2[N]; double Myu1r2[N], Myu2r2[N], Myu3r2[N], Myu4r2[N]; double Myu1i2[N], Myu2i2[N], Myu3i2[N], Myu4i2[N]; double M1r[N], M2r[N], M3r[N], M4r[N], M1a[N], M2a[N], M3a[N], M4a[N]; double M1i[N], M2i[N], M3i[N], M4i[N], M1b[N], M2b[N], M3b[N], M4b[N]; double Kr[N], Ki[N], Ka[N], Kb[N]; double Mr[N], Mi[N], Ma[N], Mb[N]; double Fr[N], Fi[N], Fa[N], Fb[N]; double sum_ka, sum_kb; double s_ka, s_kb; double ba[N], bb[N], xa[N], xb[N]; double modi_Myu2r[N], modi_Myu2i[N], modi_av[N], modi_power[N]; FILE *fp1, *fp2, *fpr, *fpi, *fkr, *fki, *fkrr, *fkii; int fft1(); int window_func(); int cmpnum(); /* check the number of input argument */ if(argc != 4) { printf("Usage: ./fft2 <No.1 input file name> <No.2 input file name> <number of frames>\n"); exit(1); } n_frame=atoi(argv[3]); // initialize arrays fp1 = fopen(argv[1], "r"); fp2 = fopen(argv[2], "r"); iter=0; flag=0; fpr = fopen("kurtosis_890_spec_ref.dat", "w"); fkr = fopen("modi_kr_890_spec_ref.txt", "w"); fki = fopen("modi_ki_890_spec_ref.txt", "w"); fkrr = fopen("kr_890_spec_ref.txt", "w"); fkii = fopen("ki_890_spec_ref.txt", "w"); while(1) { if(feof(fp1) != 0) break; else { for(i = 0; i < N; i++) { ar1[i] = 0.0; ai1[i] = 0.0; ar2[i] = 0.0; ai2[i] = 0.0; ・ ・ ・ } for(n = 0; n < n_frame; n++) { printf("\rProcessing %10d th frame",n+1); for(i = 0; i < N; i++) { fread(&ar1[i], sizeof(float), 1, fp1); //ファイル読み込み fread(&ai1[i], sizeof(float), 1, fp1); } for(i = 0; i < N; i++) { fread(&ar2[i], sizeof(float), 1, fp2); fread(&ai2[i], sizeof(float), 1, fp2); } window_func(ar1, N); //窓関数適用 window_func(ai1, N); fft1(ar1, ai1, N, iter, flag); //FFT実行 window_func(ar2, N); window_func(ai2, N); fft1(ar2, ai2, N, iter, flag); for(i = 0; i < N/2; i++) { tempr1[i] = ar1[i]; //周波数軸上での入れ替え ar1[i] = ar1[N/2+i]; ar1[N/2+i] = tempr1[i]; tempi1[i] = ai1[i]; ai1[i] = ai1[N/2+i]; ai1[N/2+i] = tempi1[i]; } for(i = 0; i < N/2; i++) { tempr2[i] = ar2[i]; ar2[i] = ar2[N/2+i]; ar2[N/2+i] = tempr2[i]; tempi2[i] = ai2[i]; ai2[i] = ai2[N/2+i]; ai2[N/2+i] = tempi2[i]; } for(i = 0; i < N; i++) { pwr1 = ar1[i]; pwi1 = ai1[i]; S1r1[i] += (double)pwr1; //和を計算 S1i1[i] += (double)pwi1; S2r1[i] += pow((double)pwr1, 2.0); //2乗和を計算 S2i1[i] += pow((double)pwi1, 2.0); S3r1[i] += pow((double)pwr1, 3.0); //3乗和を計算 S3i1[i] += pow((double)pwi1, 3.0); S4r1[i] += pow((double)pwr1, 4.0); //4乗和を計算 S4i1[i] += pow((double)pwi1, 4.0); dcnt[i] += 1.0; //繰り返した回数を計算 } for(i = 0; i < N; i++) { pwr2 = ar2[i]; pwi2 = ai2[i]; S1r2[i] += (double)pwr2; S1i2[i] += (double)pwi2; S2r2[i] += pow((double)pwr2, 2.0); S2i2[i] += pow((double)pwi2, 2.0); S3r2[i] += pow((double)pwr2, 3.0); S3i2[i] += pow((double)pwi2, 3.0); S4r2[i] += pow((double)pwr2, 4.0); S4i2[i] += pow((double)pwi2, 4.0); } } printf("\n"); for(i = 0; i < N; i++) { Myu1r1[i] = S1r1[i] / dcnt[i]; //平均を計算 Myu1i1[i] = S1i1[i] / dcnt[i]; Myu2r1[i] = S2r1[i] / dcnt[i]; Myu2i1[i] = S2i1[i] / dcnt[i]; Myu3r1[i] = S3r1[i] / dcnt[i]; Myu3i1[i] = S3i1[i] / dcnt[i]; Myu4r1[i] = S4r1[i] / dcnt[i]; Myu4i1[i] = S4i1[i] / dcnt[i]; } for(i = 0; i < N; i++) { Myu1r2[i] = S1r2[i] / dcnt[i]; Myu1i2[i] = S1i2[i] / dcnt[i]; Myu2r2[i] = S2r2[i] / dcnt[i]; Myu2i2[i] = S2i2[i] / dcnt[i]; Myu3r2[i] = S3r2[i] / dcnt[i]; Myu3i2[i] = S3i2[i] / dcnt[i]; Myu4r2[i] = S4r2[i] / dcnt[i]; Myu4i2[i] = S4i2[i] / dcnt[i]; } /* 尖度の計算 */ for(i = 0; i < N; i++) { M1r[i] = 4.0 * Myu3r1[i] * Myu1r1[i]; M1i[i] = 4.0 * Myu3i1[i] * Myu1i1[i]; M2r[i] = 6.0 * Myu2r1[i] * pow(Myu1r1[i], 2.0); M2i[i] = 6.0 * Myu2i1[i] * pow(Myu1i1[i], 2.0); M3r[i] = 3.0 * pow(Myu1r1[i], 4.0); M3i[i] = 3.0 * pow(Myu1i1[i], 4.0); M4r[i] = Myu2r1[i] - pow(Myu1r1[i], 2.0); M4i[i] = Myu2i1[i] - pow(Myu1i1[i], 2.0); M1a[i] = 4.0 * Myu3r2[i] * Myu1r2[i]; M1b[i] = 4.0 * Myu3i2[i] * Myu1i2[i]; M2a[i] = 6.0 * Myu2r2[i] * pow(Myu1r2[i], 2.0); M2b[i] = 6.0 * Myu2i2[i] * pow(Myu1i2[i], 2.0); M3a[i] = 3.0 * pow(Myu1r2[i], 4.0); M3b[i] = 3.0 * pow(Myu1i2[i], 4.0); M4a[i] = Myu2r2[i] - pow(Myu1r2[i], 2.0); M4b[i] = Myu2i2[i] - pow(Myu1i2[i], 2.0); } for(i = 0; i < N; i++) { Kr[i] = (Myu4r1[i] - M1r[i] + M2r[i] - M3r[i]) / pow(M4r[i], 2.0); Ki[i] = (Myu4i1[i] - M1i[i] + M2i[i] - M3i[i]) / pow(M4i[i], 2.0); Ka[i] = (Myu4r2[i] - M1a[i] + M2a[i] - M3a[i]) / pow(M4a[i], 2.0); Kb[i] = (Myu4i2[i] - M1b[i] + M2b[i] - M3b[i]) / pow(M4b[i], 2.0); fprintf(fkrr, "%lf ", Kr[i]); fprintf(fkii, "%lf ", Ki[i]); } fprintf(fkrr, "\n"); fprintf(fkii, "\n"); /* 尖度の計算 終了 */ /* 尖度のしきい値の計算 */ for(i = 0;i < N; i++) { Mr[i] = Kr[i] - 3.0; Mi[i] = Ki[i] - 3.0; Fr[i] = fabs(Mr[i]); Fi[i] = fabs(Mi[i]); Ma[i] = Ka[i] - 3.0; Mb[i] = Kb[i] - 3.0; Fa[i] = fabs(Ma[i]); Fb[i] = fabs(Mb[i]); } for(i = 0; i < N; i++) { sum_ka += pow(Ma[i], 2.0); sum_kb += pow(Mb[i], 2.0); } s_ka = 3.0 * sqrt(sum_ka / N); s_kb = 3.0 * sqrt(sum_kb / N); printf("s_ka=%lf, s_kb=%lf\n", s_ka, s_kb); /* 尖度のしきい値の計算 終了 */ /* RFIの判別 */ for(i = 0; i < N; i++) { if(Fr[i] > s_ka && Fi[i] > s_kb) { xa[i] = log(0.0); xb[i] = log(0.0); ba[i] = 0.0; bb[i] = 0.0; } else { xa[i] = Kr[i]; xb[i] = Ki[i]; ba[i] = Myu2r1[i]; bb[i] = Myu2i1[i]; } fprintf(fpr, "%lf ", xa[i]); fprintf(fpi, "%lf ", xb[i]); } fprintf(fpr, "\n"); fprintf(fpi, "\n"); /* RFIの判別 終了 */ /* RFI除去後のスペクトルの計算 */ for(i = 0; i < N; i++) { modi_Myu2r[i] = ba[i]; modi_Myu2i[i] = bb[i]; } for(i = 0; i < N; i++) { modi_av[i] = (modi_Myu2r[i] + modi_Myu2i[i]) / 1000000; modi_power[i] = 10 * log10(modi_av[i]); } for(i = 0; i < N; i++) { fprintf(fpr, "%lf ", modi_power[i]); } fprintf(fpr, "\n"); /* RFI除去後のスペクトルの計算 終了 */ } } fclose(fkrr); fclose(fkii); fclose(fkr); fclose(fki); fclose(fpr); fclose(fp1); fclose(fp2); return 0; } FFTと窓関数のコードは割愛
0 コメント