17 void PowerSpectrum(
float *tsIn,
float *tsOut,
int length)
20 float bprev, bnext, b0;
24 tstmp = (
float *) malloc(
sizeof(
float)*length);
25 dout = (fftwf_complex *) fftwf_malloc(
sizeof(fftwf_complex)*(length/2+1));
27 printf(
"Making fftw plan.\n");
28 p=fftwf_plan_dft_r2c_1d(length, tstmp, dout, FFTW_MEASURE);
30 printf(
"Initialize input array.\n");
31 for(i=0; i<length; i++)
36 printf(
"Executing fftw plan.\n");
39 printf(
"Calculating the absolute value.\n");
40 for(i=0; i<length/2+1; i++) {
41 tstmp[i] = cabsf(dout[i]);
44 printf(
"Interpoloating spectrum\n");
45 for(i=0; i<length/2+1; i++)
48 bnext=0.5*powf(fabsf(tstmp[i]+tstmp[i+1]), 2);
49 tsOut[i]=fmaxf(tstmp[i]*tstmp[i], bnext);
51 else if(i==length/2) {
52 bprev=0.5*powf(fabsf(tstmp[i]+tstmp[i-1]), 2);
53 tsOut[i]=fmaxf(tstmp[i]*tstmp[i], bprev);
56 bprev=0.5*powf(fabsf(tstmp[i]+tstmp[i-1]), 2);
57 bnext=0.5*powf(fabsf(tstmp[i]+tstmp[i+1]), 2);
58 b0=fmaxf(bprev, bnext);
59 tsOut[i]=fmax(b0, tstmp[i]*tstmp[i]);
63 fftwf_destroy_plan(p);
68 void delta_train(
float *tsDelta,
int length,
int period)
71 for(i=0; i<length; i++)
81 int main(
int argc,
char **argv)
86 float *tsIn = (
float *) malloc(length*
sizeof(
float));
87 float *tsOut = (
float *) malloc(length*
sizeof(
float));
90 printf(
"Making train of delta functions.\n");
91 delta_train(tsIn, length, period);
92 fout=fopen(
"input_vector.dat",
"w");
93 fwrite(tsIn,
sizeof(
float), length, fout);
96 printf(
"Calculating power spectrum.\n");
97 PowerSpectrum(tsIn, tsOut, length);
99 printf(
"Writing data to disk.\n");
101 fout=fopen(
"output_vector.dat",
"w");
105 for(i=0; i< length/2+1; ++i)
106 fprintf(fout,
"%f\n", tsOut[i]);