Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
sumhrm_v1.c
1 
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <string.h>
5 #include <assert.h>
6 
7 void main(int argc, char *argv[])
8 {
9 
10  int i=0;
11  /* the data file contains the power spectrum of dedispersed timeseries*/
12  /* The first few points (nf1 below) are set to zero as defined by */
13  /* corresponding to a period 9.9999 */
14  FILE *fp;
15  char line[256];
16 
17  float *sp; /* spectrum - a 1D array of floats*/
18  int nfolds = 5, fold; /* number of harmonic sums*/
19  int foldvals[5] = {1,2,4,8,16}; /* successive array-stretch by a factor of 2*/
20  int n,m,ii, fval;
21  float x, xdiv;
22  int lstidx;
23  int NFS;
24  int nf = 0; /* number of frequnecy channels in the spectrum*/
25  int nf1 =4; /* defined by 2*tsamp*nf/Pmax, where tsamp is the original sampling time*/
26  /* and Pmax is the largest period pulsar that we'd like to detect. */
27  /* For the test data in here, pmax=9.9999 seconds, and tsamp=54.36666 us*/
28 
29  /* The 2D array with 5 rows and nf columns. The rows contain the harmonic sums. */
30  /* first row with strech=1, or same as the original spectra, second*/
31  /* row contains 2nd harmonic sum, ie spectra streched by 2 and added to the original - */
32  /* Third row contains the 3rd harmonic sum etc. */
33 
34  /* only five harmonics - so declare the number of rows we need */
35  float **power = malloc(nfolds * sizeof(*power));
36 
37  /* following lines read the file with input spectra. Int he processes, we find out the */
38  /* number of columns and use that to dynamically allocate space*/
39  fp = fopen(argv[1], "r");
40 
41  while ( fgets(line, 256, fp) != NULL ) nf++;
42 
43  fprintf(stderr,"read %d lines of data\n",nf);
44 
45  /* we got the number of columns - declare second dimension of the array*/
46  for(fold = 0; fold < nfolds; fold++){
47  power[fold] = malloc(nf * sizeof(**power));
48  assert(power[fold] != NULL);
49  }
50 
51  /* and this array contains the original spectra*/
52  sp = (float *) malloc(nf*4);
53  assert(sp != NULL);
54 
55  /* rewind file, so that we can read from the start*/
56  fseek(fp, 0, SEEK_SET);
57 
58  i=0;
59 
60  while ( fgets(line, 250, fp) != NULL ) {
61  sp[i] = strtof(line, NULL);
62  i++;
63  };
64  fprintf(stderr,"converted %d lines of data\n",i);
65 
66  /* Harmonic sum - a direct translation of sigproc's sumhrm.f */
67  /* Fortran's array conventions are inverted as is normal in C*/
68  for (fold=0; fold< nfolds; fold++){
69  fval = foldvals[fold];
70  /* the first harmonic sum is just the copy of original array */
71  if (fval == 1) for (n=0;n<nf;n++) power[0][n] = sp[n];
72  else {
73  NFS = fval*nf1 - fval/2;
74  xdiv = 1.0/fval;
75  for (n=NFS; n< nf; n++){
76  power[fold][n] = 0.0;
77  lstidx = -1;
78  for (m=0;m<foldvals[fold];m++){
79  x = n * m * xdiv + 0.5; /* stretch by this factor before adding*/
80  ii = (int) x;
81  if ( ii > 1 && ii != lstidx) power[fold][n] = power[fold][n] + sp[ii];
82  lstidx = ii;
83  }
84  }
85  }
86  }
87 
88  /* write out the first 20000 values*/
89  for (ii=0;ii<20000;ii++){
90  fprintf(stdout, "%07d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
91  ii, sp[ii],power[0][ii],power[1][ii],power[2][ii],power[3][ii],power[4][ii]);
92  printf("\n");
93  }
94 
95  /* release the memory area and close data file*/
96  for(fold = 0; fold < nfolds; fold++) free(power[fold]);
97 
98  free(power);
99 
100  fclose(fp);
101 }