Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
TDT_Powerspectrum.c
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <complex.h>
5 #include <fftw3.h>
6 
7 /*cmodel for the TDT PWFT block.
8 This function makes use of FFTW.
9 Note that because I included complex.h, the output data time from the
10 FFTW plan/execute is a C complex time, not a double[2].
11 I also used the float versions of the functions,
12 which must be linked with -lfftw3f (not -lfftw3)
13 e.g.:
14 gcc TDT_Powerspectrum.c -o powerspectrum -lfftw3f -lm
15 */
16 
17 void PowerSpectrum(float *tsIn, float *tsOut, int length)
18 {
19  int i;
20  float bprev, bnext, b0;
21  fftwf_plan p;
22  fftwf_complex *dout;
23  float *tstmp;
24  tstmp = (float *) malloc(sizeof(float)*length);
25  dout = (fftwf_complex *) fftwf_malloc(sizeof(fftwf_complex)*(length/2+1));
26 
27  printf("Making fftw plan.\n");
28  p=fftwf_plan_dft_r2c_1d(length, tstmp, dout, FFTW_MEASURE);
29 
30  printf("Initialize input array.\n");
31  for(i=0; i<length; i++)
32  {
33  tstmp[i]=tsIn[i];
34  }
35 
36  printf("Executing fftw plan.\n");
37  fftwf_execute(p);
38 
39  printf("Calculating the absolute value.\n");
40  for(i=0; i<length/2+1; i++) {
41  tstmp[i] = cabsf(dout[i]); //NOTE: cabsf defined in complex.h
42  }
43 
44  printf("Interpoloating spectrum\n");
45  for(i=0; i<length/2+1; i++)
46  {
47  if(i==0) {
48  bnext=0.5*powf(fabsf(tstmp[i]+tstmp[i+1]), 2);
49  tsOut[i]=fmaxf(tstmp[i]*tstmp[i], bnext);
50  }
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);
54  }
55  else {
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]);
60  }
61  }
62 
63  fftwf_destroy_plan(p);
64  fftwf_free(dout);
65  free(tstmp);
66 }
67 
68 void delta_train(float *tsDelta, int length, int period)
69 {
70  int i;
71  for(i=0; i<length; i++)
72  {
73  if(i%period==0) {
74  tsDelta[i]=1.0;
75  } else {
76  tsDelta[i]=0.0;
77  }
78  }
79 }
80 
81 int main(int argc, char **argv)
82 {
83  int length=32768;
84  int period=256; //In samples
85 
86  float *tsIn = (float *) malloc(length*sizeof(float));
87  float *tsOut = (float *) malloc(length*sizeof(float));
88  FILE *fout;
89 
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);
94  fclose(fout);
95 
96  printf("Calculating power spectrum.\n");
97  PowerSpectrum(tsIn, tsOut, length);
98 
99  printf("Writing data to disk.\n");
100 
101  fout=fopen("output_vector.dat", "w");
102  //fwrite(tsOut, sizeof(float), length/2+1, fout);
103 
104  int i;
105  for(i=0; i< length/2+1; ++i)
106  fprintf(fout, "%f\n", tsOut[i]);
107 
108  fclose(fout);
109 
110  free(tsIn);
111  free(tsOut);
112 
113 }
114