text: /*
   * @file yapp_ft.c
   * Program to compute [PFB +] FFT on baseband (raw voltage time series) data
   *
   * @verbatim
   * Usage: yapp_viewdata [options] <data-file>
   *     -h  --help                           Display this usage information
   *     -s  --skip <time>                    The length of data in seconds, to be
   *                                          skipped
   *                                          (default is 0 s)
   *     -p  --proc <time>                    The length of data in seconds, to be
   *                                          processed
   *                                          (default is all)
   *     -l  --tsamp <tsamp>                  Sampling time in s
   *     -f  --centre-freq <freq>             Centre frequency of observing band
   *     -b  --pfb                            Do polyphase filter bank
   *     -t  --ntaps <ntaps>                  Number of taps in the PFB
   *                                          (default is 8)
   *     -n  --nfft <nfft>                    Number of points in FFT
   *                                          (default is 1024)
   *     -a  --nacc <acc-len>                 Number of spectra to accumulate
   *                                          (default is 1)
   *     -o  --obs-site <site-name>           Observatory name
   *     -r  --pulsar <psr-name>              Pulsar name
   *     -g  --graphics                       Turn on plotting
   *     -i  --invert                         Invert the background and foreground
   *                                          colours in plots
   *     -e  --non-interactive                Run in non-interactive mode
   *     -v  --version                        Display the version @endverbatim
   *
   * @author Jayanth Chennamangalam
   * @date 2012.07.18
   */
  
  #include "yapp.h"
  #include "yapp_ft.h"
  
  /**
   * The build version string, maintained in the file version.c, which is
   * generated by makever.c.
   */
  extern const char *g_pcVersion;
  
  /* PGPLOT device ID */
  extern int g_iPGDev;
  
  /* the following are global only to enable cleaning up in case of abnormal
     termination, such as those triggered by SIGINT or SIGTERM */
  float *g_pfXAxis = NULL;
  PFB_DATA g_astPFBData[NUM_TAPS] = {{0}};
  int g_iPFBReadIdx = 0;
  int g_iPFBWriteIdx = 0;
  fftwf_complex *g_pfcFFTInX = NULL;
  fftwf_complex *g_pfcFFTOutX = NULL;
  fftwf_plan g_stPlanX = {0};
  fftwf_complex *g_pfcFFTInY = NULL;
  fftwf_complex *g_pfcFFTOutY = NULL;
  fftwf_plan g_stPlanY = {0};
  float *g_pfSumPow = NULL;
  float *g_pfSumPowX = NULL;
  float *g_pfSumPowY = NULL;
  float *g_pfSumStokesRe = NULL;
  float *g_pfSumStokesIm = NULL;
  char g_cDoPFB = YAPP_FALSE;
  char g_acFileData[LEN_GENSTRING] = {0};
  char g_acFileCoeff[LEN_GENSTRING] = {0};
  int g_iFileData = 0;
  float* g_pfPFBCoeff = NULL;
  
  int main(int argc, char *argv[])
  {
      double dDataSkipTime = 0.0;
      double dDataProcTime = 0.0;
      YUM_t stYUM = {{0}};
      double dTSampInSec = 0.0;   /* holds sampling time in s */
      int iBlockSize = 0;
      int iNumAcc = DEF_ACC;
      int iFileSpec = 0;
      int iFileSpecHdr = 0;
      int iSpecCount = 0;
      struct stat stFileStats = {0};
      long lBytesToSkip = 0;
      long lBytesToProc = 0;
      int iTimeSampsToSkip = 0;
      int iTimeSampsToProc = 0;
      char cIsFil = YAPP_TRUE;
      int iNTaps = 1;                       /* 1 if no PFB, NUM_TAPS if PFB */
      int iNFFT = DEF_NFFT;
      int iNumReads = 0;
      int iReadBlockCount = 0;
      float fFreqCen = 0.0;
      char cIsFirst = YAPP_TRUE;
      char cIsLastBlock = YAPP_FALSE;
      int iRet = YAPP_RET_SUCCESS;
      float fDataMin = 0.0;
      float fDataMax = 0.0;
      float fColMin = 0.0;
      float fColMax = 0.0;
      float fButX = 0.0;
      float fButY = 0.0;
      char cCurChar = 0;
      long int lDataSizeTotal = 0;
      float fSampSize = 0.0; 
      int iTimeSamps = 0;
      char acFileSpec[LEN_GENSTRING] = {0};
      char acFileSpecHdr[LEN_GENSTRING] = {0};
      char acHdrBuf[LEN_GENSTRING] = {0};
      char acSite[LEN_GENSTRING] = {0};
      char acPulsar[MAX_LEN_PSRNAME] = {0};
      int i = 0;
      int j = 0;
      char cHasGraphics = YAPP_FALSE;
      int iInvCols = YAPP_FALSE;
      char cIsNonInteractive = YAPP_FALSE;
      const char *pcProgName = NULL;
      int iNextOpt = 0;
      /* valid short options */
      const char* const pcOptsShort = "hs:p:l:f:bt:n:a:o:r:giev";
      /* valid long options */
      const struct option stOptsLong[] = {
          { "help",                   0, NULL, 'h' },
          { "skip",                   1, NULL, 's' },
          { "proc",                   1, NULL, 'p' },
          { "tsamp",                  1, NULL, 'l' },
          { "centre-freq",            1, NULL, 'f' },
          { "pfb",                    0, NULL, 'b' },
          { "ntaps",                  1, NULL, 't' },
          { "nfft",                   1, NULL, 'n' },
          { "nacc",                   1, NULL, 'a' },
          { "obs-site",               1, NULL, 'o' }, 
          { "pulsar",                 1, NULL, 'r' },
          { "graphics",               0, NULL, 'g' },
          { "invert",                 0, NULL, 'i' },
          { "non-interactive",        0, NULL, 'e' },
          { "version",                0, NULL, 'v' },
          { NULL,                     0, NULL, 0   }
      };
  
      /* get the filename of the program from the argument list */
      pcProgName = argv[0];
  
      /* parse the input */
      do
      {
          iNextOpt = getopt_long(argc, argv, pcOptsShort, stOptsLong, NULL);
          switch (iNextOpt)
          {
              case 'h':   /* -h or --help */
                  /* print usage info and terminate */
                  PrintUsage(pcProgName);
                  return YAPP_RET_SUCCESS;
  
              case 's':   /* -s or --skip */
                  /* set option */
                  dDataSkipTime = atof(optarg);
                  break;
  
              case 'p':   /* -p or --proc */
                  /* set option */
                  dDataProcTime = atof(optarg);
                  break;
  
              case 'l':   /* -l or --tsamp */
                  /* set option */
                  dTSampInSec = atof(optarg);
                  break;
  
              case 'f':   /* -f or --centre-freq */
                  /* set option */
                  fFreqCen = atof(optarg);
                  break;
  
              case 'b':   /* -b or --pfb */
                  /* set option */
                  g_cDoPFB = YAPP_TRUE;
                  break;
  
              case 't':   /* -t or --ntaps */
                  /* set option */
                  iNTaps = atoi(optarg);
                  break;
  
              case 'n':   /* -n or --nfft */
                  /* set option */
                  iNFFT = atoi(optarg);
                  break;
  
              case 'a':   /* -a or --nacc */
                  /* set option */
                  iNumAcc = (int) atoi(optarg);
                  break;
  
              case 'o':   /* -o or --obs-site */
                  /* set option */
                  (void) strncpy(acSite, optarg, LEN_GENSTRING);
                  break;
  
              case 'r':   /* -r or --pulsar */
                  /* set option */
                  (void) strncpy(acPulsar, optarg, MAX_LEN_PSRNAME);
                  break;
  
              case 'g':   /* -g or --graphics */
                  /* set option */
                  cHasGraphics = YAPP_TRUE;
                  break;
  
              case 'i':  /* -i or --invert */
                  /* set option */
                  iInvCols = YAPP_TRUE;
                  break;
  
              case 'e':  /* -e or --non-interactive */
                  /* set option */
                  cIsNonInteractive = YAPP_TRUE;
                  break;
  
              case 'v':   /* -v or --version */
                  /* display the version */
                  (void) printf("%s\n", g_pcVersion);
                  return YAPP_RET_SUCCESS;
  
              case '?':   /* user specified an invalid option */
                  /* print usage info and terminate with error */
                  (void) fprintf(stderr, "ERROR: Invalid option!\n");
                  PrintUsage(pcProgName);
                  return YAPP_RET_ERROR;
  
              case -1:    /* done with options */
                  break;
  
              default:    /* unexpected */
                  assert(0);
          }
      } while (iNextOpt != -1);
  
      /* no arguments */
      if (argc <= optind)
      {
          (void) fprintf(stderr, "ERROR: Input file not specified!\n");
          PrintUsage(pcProgName);
          return YAPP_RET_ERROR;
      }
  
      if (0.0 == dTSampInSec)
      {
          (void) fprintf(stderr, "ERROR: Sampling interval not specified!\n");
          PrintUsage(pcProgName);
          return YAPP_RET_ERROR;
      }
  
      if (0.0 == fFreqCen)
      {
          (void) fprintf(stderr, "ERROR: Centre frequency not specified!\n");
          PrintUsage(pcProgName);
          return YAPP_RET_ERROR;
      }
  
      if ((g_cDoPFB) && (1 == iNTaps))
      {
          /* set default number of taps */
          iNTaps = NUM_TAPS;
      }
  
      /* set the read block size */
      iBlockSize = iNFFT;
  
      /* register the signal-handling function */
      iRet = YAPP_RegisterSignalHandlers();
      if (iRet != YAPP_RET_SUCCESS)
      {
          (void) fprintf(stderr,
                         "ERROR: Handler registration failed!\n");
          return YAPP_RET_ERROR;
      }
  
      /* get the input filename */
      (void) strncpy(g_acFileData, argv[optind], LEN_GENSTRING); 
  
  
      iRet = stat(g_acFileData, &stFileStats);
      if (iRet != YAPP_RET_SUCCESS)
      {
          (void) fprintf(stderr,
                         "ERROR: Failed to stat %s: %s!\n",
                         g_acFileData,
                         strerror(errno));
          CleanUp(iNTaps);
          return YAPP_RET_ERROR;
      }
      lDataSizeTotal = (long) stFileStats.st_size;
      fSampSize = NUM_BYTES_PER_SAMP * sizeof(char); 
      iTimeSamps = (int) floor(((double) lDataSizeTotal) / fSampSize);
  
      /* calculate bytes to skip and read */
      if (0.0 == dDataProcTime)
      {
          dDataProcTime = (iTimeSamps * dTSampInSec) - dDataSkipTime;
      }
      /* check if the input time duration is less than the length of the
         data */
      else if (dDataProcTime > (iTimeSamps * dTSampInSec))
      {
          (void) fprintf(stderr,
                         "WARNING: Input time is longer than length of "
                         "data!\n");
      }
  
      lBytesToSkip = (long) floor(dDataSkipTime / dTSampInSec)
                                                      /* number of samples */
                            * fSampSize;
      if (fSampSize >= 2.0)
      {
          /* will be true for fSampSize = 0.5 and 1.0 */
          assert(0 == lBytesToSkip % (int) fSampSize);
      }
      lBytesToProc = (long) floor(dDataProcTime / dTSampInSec)
                                                      /* number of samples */
                            * fSampSize;
      if (fSampSize >= 2.0)
      {
          /* will be true for fSampSize = 0.5 and 1.0 */
          assert(0 == lBytesToProc % (int) fSampSize);
      }
  
      if (lBytesToSkip >= stYUM.lDataSizeTotal)
      {
          (void) fprintf(stderr,
                         "ERROR: Data to be skipped is greater than or equal to "
                         "the size of the file!\n");
          YAPP_CleanUp();
          return YAPP_RET_ERROR;
      }
  
      if ((lBytesToSkip + lBytesToProc) > lDataSizeTotal)
      {
          (void) printf("WARNING: Total data to be read (skipped and processed) "
                        "is more than the size of the file! ");
          lBytesToProc = lDataSizeTotal - lBytesToSkip;
          (void) printf("Newly calculated size of data to be processed: %ld "
                        "bytes\n",
                        lBytesToProc);
      }
  
      iTimeSampsToSkip = (int) (lBytesToSkip / fSampSize);
      (void) printf("Skipping\n"
                    "    %ld of %ld bytes\n"
                    "    %d of %d time samples\n"
                    "    %.10g of %.10g seconds\n",
                    lBytesToSkip,
                    lDataSizeTotal,
                    iTimeSampsToSkip,
                    iTimeSamps,
                    (iTimeSampsToSkip * dTSampInSec),
                    (iTimeSamps * dTSampInSec));
  
      /* truncate lBytesToProc to (X * iNTaps * iNFFT * NUM_BYTES_PER_SAMP)
         bytes where X is the largest possible integer, if more */
      float fBlocks = (float) lBytesToProc
                      / (iNTaps * iNFFT * NUM_BYTES_PER_SAMP);
      int iBlocks = (int) fBlocks;
      if (fBlocks != (float) iBlocks)
      {
          size_t iDiffBytes = (size_t) ((fBlocks - (float) iBlocks)
                                        * iNTaps * iNFFT * NUM_BYTES_PER_SAMP);
          lBytesToProc -= iDiffBytes;
      }
  
      iTimeSampsToProc = (int) (lBytesToProc / fSampSize);
      iNumReads = (int) floorf(((float) iTimeSampsToProc) / iBlockSize);
  
      /* initialise */
      iRet = InitPFB(iNTaps, iNFFT);
      if (iRet != YAPP_RET_SUCCESS)
      {
          (void) fprintf(stderr, "ERROR! InitPFB failed!\n");
          CleanUp(iNTaps);
          return YAPP_RET_ERROR;
      }
  
      g_iFileData = open(g_acFileData, O_RDONLY);
      if (g_iFileData < YAPP_RET_SUCCESS)
      {
          (void) fprintf(stderr,
                         "ERROR! Opening data file %s failed! %s.\n",
                         g_acFileData,
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
  
      (void) printf("Processing\n"
                    "    %ld of %ld bytes\n"
                    "    %d of %d time samples\n"
                    "    %.10g of %.10g seconds\n"
                    "in %d reads with block size %d time samples...\n",
                    lBytesToProc,
                    lDataSizeTotal,
                    iTimeSampsToProc,
                    iTimeSamps,
                    (iTimeSampsToProc * dTSampInSec),
                    (iTimeSamps * dTSampInSec),
                    iNumReads,
                    iBlockSize);
  
      /* open output file */
      if (cIsFil)
      {
          /* create output file */
          /* this will be freed by the YAPP garbage collector */
          char* pcFilename = YAPP_GetFilenameWithExtFromPath(g_acFileData);
          (void) strcpy(acFileSpec, pcFilename);
          (void) strcat(acFileSpec, EXT_FIL);
          iFileSpec = open(acFileSpec,
                           O_CREAT | O_TRUNC | O_WRONLY,
                           S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);
          if (iFileSpec < YAPP_RET_SUCCESS)
          {
              (void) fprintf(stderr, "ERROR: Opening spectrum file failed!\n");
              CleanUp(iNTaps);
              return YAPP_RET_ERROR;
          }
      
          /* open fil header file */
          /* this will be freed by the YAPP garbage collector */
          pcFilename = YAPP_GetFilenameWithExtFromPath(g_acFileData);
          (void) strcpy(acFileSpecHdr, pcFilename);
          (void) strcat(acFileSpecHdr, EXT_YM);
          iFileSpecHdr = open(acFileSpecHdr,
                              O_CREAT | O_TRUNC | O_WRONLY,
                              S_IRUSR | S_IWUSR | S_IRGRP | S_IROTH);
          if (iFileSpecHdr < YAPP_RET_SUCCESS)
          {
              (void) fprintf(stderr,
                             "ERROR: Opening spectrum header file failed!\n");
              CleanUp(iNTaps);
              return YAPP_RET_ERROR;
          }
  
          /* populate the output metadata structure */
          (void) strncpy(stYUM.acSite, acSite, LEN_GENSTRING);
          (void) strncpy(stYUM.acPulsar, acPulsar, MAX_LEN_PSRNAME);
          stYUM.dTSamp = dTSampInSec * 1e3 * iNFFT * iNumAcc;
          stYUM.iTimeSamps = (int) ((float) iTimeSamps / (iNFFT * iNumAcc));
          stYUM.iNumChans = iNFFT;
          stYUM.iNumGoodChans = iNFFT;
          stYUM.fChanBW = (((float) 1) / (dTSampInSec * iNFFT)) / 1e6;
          stYUM.fFMin = fFreqCen - ((stYUM.iNumChans / 2) * stYUM.fChanBW);
          stYUM.fFMax = fFreqCen + ((stYUM.iNumChans / 2) * stYUM.fChanBW);
          stYUM.iNumBands = 1;
          stYUM.iNumBadTimes = 0;
          stYUM.iNumBits = 32;
          stYUM.iNumIFs = 1;
  
          /* write header info */
          (void) sprintf(acHdrBuf,
                         "Observing site                    : %s\n",
                         stYUM.acSite);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Field name                        : %s\n",
                         stYUM.acPulsar);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Sampling interval                 : %.10g ms\n",
                         stYUM.dTSamp);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Number of channels                : %d\n",
                         stYUM.iNumChans);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Number of good channels           : %d\n",
                         stYUM.iNumGoodChans);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Channel bandwidth                 : %.10g MHz\n",
                         stYUM.fChanBW);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Lowest frequency                  : %.10g MHz\n",
                         stYUM.fFMin);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Highest frequency                 : %.10g MHz\n",
                         stYUM.fFMax);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Estimated number of bands         : %d\n",
                         stYUM.iNumBands);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Number of bad time sections       : %d\n",
                         stYUM.iNumBadTimes);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Number of bits per sample         : %d\n",
                         stYUM.iNumBits);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "Number of IFs                     : %d\n",
                         stYUM.iNumIFs);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf, "Duration of data in\n");
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "    Bytes                         : %ld\n",
                         (lDataSizeTotal / (iNFFT * iNumAcc)));
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "    Time samples                  : %d\n",
                         stYUM.iTimeSamps);
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
          (void) sprintf(acHdrBuf,
                         "    Time                          : %g s\n",
                         (stYUM.iTimeSamps * (stYUM.dTSamp * 1e-3)));
          (void) write(iFileSpecHdr, acHdrBuf, strlen(acHdrBuf));
  
          (void) close(iFileSpecHdr);
      }
  
      if (1 == iNumReads)
      {
          cIsLastBlock = YAPP_TRUE;
      }
  
      /* skip data, if any are to be skipped */
      (void) lseek(g_iFileData, lBytesToSkip, SEEK_SET);
  
      if (cHasGraphics)
      {
          /* open the PGPLOT graphics device */
          g_iPGDev = cpgopen(PG_DEV);
          if (g_iPGDev <= 0)
          {
              (void) fprintf(stderr,
                             "ERROR: Opening graphics device %s failed!\n",
                             PG_DEV);
              CleanUp(iNTaps);
              return YAPP_RET_ERROR;
          }
  
          /* set the background colour to white and the foreground colour to
             black, if user requires so */
          if (YAPP_TRUE == iInvCols)
          {
              cpgscr(0, 1.0, 1.0, 1.0);
              cpgscr(1, 0.0, 0.0, 0.0);
          }
  
          /* set up the plot's X-axis */
          g_pfXAxis = (float *) YAPP_Malloc(iBlockSize,
                                            sizeof(float),
                                            YAPP_FALSE);
          if (NULL == g_pfXAxis)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation for X-axis failed! %s!\n",
                             strerror(errno));
              CleanUp(iNTaps);
              return YAPP_RET_ERROR;
          }
          /* load the X-axis frequency axis */
          for (i = 0; i < iBlockSize; ++i)
          {
              /* frequency axis in MHz */
              g_pfXAxis[i] = ((float) i * 1e-6) / (dTSampInSec * iNFFT);
          }
      }
  
      while (iNumReads > 0)
      {
          /* read data */
          (void) printf("\rReading data block %d.", iReadBlockCount);
          (void) fflush(stdout);
          /* read data from input buffer */
          iRet = ReadData(cIsFirst, iNTaps, iNFFT);
          if (YAPP_RET_ERROR == iRet)
          {
              (void) fprintf(stderr, "ERROR: Reading data failed!\n");
              CleanUp(iNTaps);
              return YAPP_RET_ERROR;
          }
          else if (YAPP_RET_READDONE == iRet)
          {
              break;
          }
          cIsFirst = YAPP_FALSE;
          --iNumReads;
          ++iReadBlockCount;
  
          if (g_cDoPFB)
          {
              /* do pfb */
              (void) DoPFB(iNTaps, iNFFT);
          }
          else
          {
              /* copy data for FFT */
              (void) CopyDataForFFT(iNFFT);
          }
  
          /* do fft */
          iRet = DoFFT();
          if (iRet != YAPP_RET_SUCCESS)
          {
              (void) fprintf(stderr, "ERROR! FFT failed!\n");
              (void) close(iFileSpec);
              CleanUp(iNTaps);
              return YAPP_RET_ERROR;
          }
  
          if (cIsFil)
          {
              for (i = 0; i < iNFFT; ++i)
              {
                  g_pfSumPow[i] += (g_pfcFFTOutX[i][0] * g_pfcFFTOutX[i][0])
                                   + (g_pfcFFTOutX[i][1] * g_pfcFFTOutX[i][1])
                                   + (g_pfcFFTOutY[i][0] * g_pfcFFTOutY[i][0])
                                   + (g_pfcFFTOutY[i][1] * g_pfcFFTOutY[i][1]);
              }
          }
          else
          {
              /* accumulate power x, power y, stokes */
              for (i = 0; i < iNFFT; ++i)
              {
                  /* Re(X)^2 + Im(X)^2 */
                  g_pfSumPowX[i] += (g_pfcFFTOutX[i][0] * g_pfcFFTOutX[i][0])
                                    + (g_pfcFFTOutX[i][1] * g_pfcFFTOutX[i][1]);
                  /* Re(Y)^2 + Im(Y)^2 */
                  g_pfSumPowY[i] += (g_pfcFFTOutY[i][0] * g_pfcFFTOutY[i][0])
                                    + (g_pfcFFTOutY[i][1] * g_pfcFFTOutY[i][1]);
                  /* Re(XY*) */
                  g_pfSumStokesRe[i] += (g_pfcFFTOutX[i][0] * g_pfcFFTOutY[i][0])
                      + (g_pfcFFTOutX[i][1] * g_pfcFFTOutY[i][1]);
                  /* Im(XY*) */
                  g_pfSumStokesIm[i] += (g_pfcFFTOutX[i][1] * g_pfcFFTOutY[i][0])
                      - (g_pfcFFTOutX[i][0] * g_pfcFFTOutY[i][1]);
              }
          }
          ++iSpecCount;
          if (iSpecCount == iNumAcc)
          {
              if (cIsFil)
              {
                  (void) write(iFileSpec, g_pfSumPow, iNFFT * sizeof(float));
  
                  if (cHasGraphics)
                  {
                      fDataMin = g_pfSumPow[0];
                      fDataMax = g_pfSumPow[0];
                      for (j = 0; j < iBlockSize; ++j)
                      {
                          if (g_pfSumPow[j] < fDataMin)
                          {
                              fDataMin = g_pfSumPow[j];
                          }
                          if (g_pfSumPow[j] > fDataMax)
                          {
                              fDataMax = g_pfSumPow[j];
                          }
                      }
  
                      #ifdef DEBUG
                      (void) printf("Minimum value of data             : %g\n",
                                    fDataMin);
                      (void) printf("Maximum value of data             : %g\n",
                                    fDataMax);
                      #endif
  
                      fColMin = fDataMin;
                      fColMax = fDataMax;
  
                      if (!(cIsLastBlock))
                      {
                          /* erase just before plotting, to reduce flicker */
                          cpgeras();
                      }
  
                      cpgsvp(PG_VP_ML, PG_VP_MR, PG_VP_MB, PG_VP_MT);
                      cpgswin(g_pfXAxis[0],
                              g_pfXAxis[iBlockSize-1],
                              fColMin,
                              fColMax);
                      cpglab("Frequency (MHz)", "Total Power", "Power Spectrum");
                      cpgbox("BCNST", 0.0, 0, "BCNST", 0.0, 0);
                      cpgsci(PG_CI_PLOT);
                      cpgline(iNFFT, g_pfXAxis, g_pfSumPow);
                      cpgsci(PG_CI_DEF);
  
                      if (!(cIsLastBlock))
                      {
                          if (!(cIsNonInteractive))
                          {
                              /* draw the 'next' and 'exit' buttons */
                              cpgsvp(PG_VP_BUT_ML, PG_VP_BUT_MR, PG_VP_BUT_MB, PG_VP_BUT_MT);
                              cpgswin(PG_BUT_L, PG_BUT_R, PG_BUT_B, PG_BUT_T);
                              cpgsci(PG_BUT_FILLCOL); /* set the fill colour */
                              cpgrect(PG_BUTNEXT_L, PG_BUTNEXT_R, PG_BUTNEXT_B, PG_BUTNEXT_T);
                              cpgrect(PG_BUTEXIT_L, PG_BUTEXIT_R, PG_BUTEXIT_B, PG_BUTEXIT_T);
                              cpgsci(0);  /* set colour index to white */
                              cpgtext(PG_BUTNEXT_TEXT_L, PG_BUTNEXT_TEXT_B, "Next");
                              cpgtext(PG_BUTEXIT_TEXT_L, PG_BUTEXIT_TEXT_B, "Exit");
  
                              fButX = (PG_BUTNEXT_R - PG_BUTNEXT_L) / 2;
                              fButY = (PG_BUTNEXT_T - PG_BUTNEXT_B) / 2;
  
                              while (YAPP_TRUE)
                              {
                                  iRet = cpgcurs(&fButX, &fButY, &cCurChar);
                                  if (0 == iRet)
                                  {
                                      (void) fprintf(stderr,
                                                     "WARNING: "
                                                     "Reading cursor parameters failed!\n");
                                      break;
                                  }
  
                                  if (((fButX >= PG_BUTNEXT_L) && (fButX <= PG_BUTNEXT_R))
                                      && ((fButY >= PG_BUTNEXT_B) && (fButY <= PG_BUTNEXT_T)))
                                  {
                                      /* animate button click */
                                      cpgsci(PG_BUT_FILLCOL);
                                      cpgtext(PG_BUTNEXT_TEXT_L, PG_BUTNEXT_TEXT_B, "Next");
                                      cpgsci(0);  /* set colour index to white */
                                      cpgtext(PG_BUTNEXT_CL_TEXT_L, PG_BUTNEXT_CL_TEXT_B, "Next");
                                      (void) usleep(PG_BUT_CL_SLEEP);
                                      cpgsci(PG_BUT_FILLCOL); /* set colour index to fill
                                                                 colour */
                                      cpgtext(PG_BUTNEXT_CL_TEXT_L, PG_BUTNEXT_CL_TEXT_B, "Next");
                                      cpgsci(0);  /* set colour index to white */
                                      cpgtext(PG_BUTNEXT_TEXT_L, PG_BUTNEXT_TEXT_B, "Next");
                                      cpgsci(1);  /* reset colour index to black */
                                      (void) usleep(PG_BUT_CL_SLEEP);
  
                                      break;
                                  }
                                  else if (((fButX >= PG_BUTEXIT_L) && (fButX <= PG_BUTEXIT_R))
                                      && ((fButY >= PG_BUTEXIT_B) && (fButY <= PG_BUTEXIT_T)))
                                  {
                                      /* animate button click */
                                      cpgsci(PG_BUT_FILLCOL);
                                      cpgtext(PG_BUTEXIT_TEXT_L, PG_BUTEXIT_TEXT_B, "Exit");
                                      cpgsci(0);  /* set colour index to white */
                                      cpgtext(PG_BUTEXIT_CL_TEXT_L, PG_BUTEXIT_CL_TEXT_B, "Exit");
                                      (void) usleep(PG_BUT_CL_SLEEP);
                                      cpgsci(PG_BUT_FILLCOL); /* set colour index to fill
                                                                 colour */
                                      cpgtext(PG_BUTEXIT_CL_TEXT_L, PG_BUTEXIT_CL_TEXT_B, "Exit");
                                      cpgsci(0);  /* set colour index to white */
                                      cpgtext(PG_BUTEXIT_TEXT_L, PG_BUTEXIT_TEXT_B, "Exit");
                                      cpgsci(1);  /* reset colour index to black */
                                      (void) usleep(PG_BUT_CL_SLEEP);
  
                                      (void) close(iFileSpec);
                                      (void) printf("\nDONE!\n");
                                      cpgclos();
                                      CleanUp(iNTaps);
                                      return YAPP_RET_SUCCESS;
                                  }
                              }
                          }
                          else
                          {
                              /* pause before erasing */
                              (void) usleep(PG_PLOT_SLEEP);
                          }
                      }
                  }
              }
              else
              {
                  (void) write(iFileSpec, g_pfSumPowX, iNFFT * sizeof(float));
                  (void) write(iFileSpec, g_pfSumPowY, iNFFT * sizeof(float));
                  (void) write(iFileSpec, g_pfSumStokesRe, iNFFT * sizeof(float));
                  (void) write(iFileSpec, g_pfSumStokesIm, iNFFT * sizeof(float));
              }
  
              /* reset time */
              iSpecCount = 0;
              /* zero accumulators */
              if (cIsFil)
              {
                  (void) memset(g_pfSumPow, '\0', iNFFT * sizeof(float));
              }
              else
              {
                  (void) memset(g_pfSumPowX, '\0', iNFFT * sizeof(float));
                  (void) memset(g_pfSumPowY, '\0', iNFFT * sizeof(float));
                  (void) memset(g_pfSumStokesRe, '\0', iNFFT * sizeof(float));
                  (void) memset(g_pfSumStokesIm, '\0', iNFFT * sizeof(float));
              }
          }
  
          if (1 == iNumReads)
          {
              cIsLastBlock = YAPP_TRUE;
          }
      }
  
      (void) close(iFileSpec);
      (void) printf("\nDONE!\n");
      if (cHasGraphics)
      {
          cpgclos();
      }
      CleanUp(iNTaps);
  
      return YAPP_RET_SUCCESS;
  }
  
  /* function that creates the FFT plan, allocates memory, initialises counters,
     etc. */
  int InitPFB(int iNTaps, int iNFFT)
  {
      int i = 0;
      int iRet = YAPP_RET_SUCCESS;
  
      if (g_cDoPFB)
      {
          int iFileCoeff = 0;
  
          g_pfPFBCoeff = (float *) YAPP_Malloc((size_t) iNTaps * iNFFT,
                                               sizeof(float),
                                               YAPP_FALSE);
          if (NULL == g_pfPFBCoeff)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s.\n",
                             strerror(errno));
              return YAPP_RET_ERROR;
          }
  
          /* read filter coefficients */
          /* build file name */
          (void) sprintf(g_acFileCoeff,
                         "%s_%d_%d_%d%s",
                         FILE_COEFF_PREFIX,
                         iNTaps,
                         iNFFT,
                         DEF_NUM_SUBBANDS,
                         FILE_COEFF_SUFFIX);
          iFileCoeff = open(g_acFileCoeff, O_RDONLY);
          if (iFileCoeff < YAPP_RET_SUCCESS)
          {
              (void) fprintf(stderr,
                             "ERROR: Opening filter coefficients file %s "
                             "failed! %s.\n",
                             g_acFileCoeff,
                             strerror(errno));
              return YAPP_RET_ERROR;
          }
  
          iRet = read(iFileCoeff,
                      g_pfPFBCoeff,
                      iNTaps * iNFFT * sizeof(float));
          if (iRet != (iNTaps * iNFFT * sizeof(float)))
          {
              (void) fprintf(stderr,
                             "ERROR: Reading filter coefficients failed! %s.\n",
                             strerror(errno));
              (void) close(iFileCoeff);
              return YAPP_RET_ERROR;
          }
          (void) close(iFileCoeff);
      }
  
      /* allocate memory for data array contents */
      for (i = 0; i < iNTaps; ++i)
      {
          g_astPFBData[i].pcData = (signed char*) YAPP_Malloc((size_t) iNFFT * NUM_BYTES_PER_SAMP,
                                                              sizeof(signed char),
                                                              YAPP_FALSE);
          if (NULL == g_astPFBData[i].pcData)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s.\n",
                             strerror(errno));
              return YAPP_RET_ERROR;
          }
          g_astPFBData[i].pfcDataX = (fftwf_complex *) fftwf_malloc(iNFFT
                                                        * sizeof(fftwf_complex));
          if (NULL == g_astPFBData[i].pfcDataX)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s.\n",
                             strerror(errno));
              return YAPP_RET_ERROR;
          }
          g_astPFBData[i].pfcDataY = (fftwf_complex *) fftwf_malloc(iNFFT
                                                        * sizeof(fftwf_complex));
          if (NULL == g_astPFBData[i].pfcDataY)
          {
              (void) fprintf(stderr,
                             "ERROR: Memory allocation failed! %s.\n",
                             strerror(errno));
              return YAPP_RET_ERROR;
          }
          if (i != (iNTaps - 1))
          {
              g_astPFBData[i].iNextIdx = i + 1;
          }
          else
          {
              g_astPFBData[i].iNextIdx = 0;
          }
      }
  
      g_pfcFFTInX = (fftwf_complex *) fftwf_malloc(iNFFT
                                                   * sizeof(fftwf_complex));
      if (NULL == g_pfcFFTInX)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
      g_pfcFFTInY = (fftwf_complex *) fftwf_malloc(iNFFT
                                                   * sizeof(fftwf_complex));
      if (NULL == g_pfcFFTInY)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
      g_pfcFFTOutX = (fftwf_complex *) fftwf_malloc(iNFFT
                                                    * sizeof(fftwf_complex));
      if (NULL == g_pfcFFTOutX)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
      g_pfcFFTOutY = (fftwf_complex *) fftwf_malloc(iNFFT
                                                    * sizeof(fftwf_complex));
      if (NULL == g_pfcFFTOutY)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
  
      g_pfSumPow = (float *) YAPP_Malloc((size_t) iNFFT,
                                         sizeof(float),
                                         YAPP_TRUE);
      if (NULL == g_pfSumPow)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
      g_pfSumPowX = (float *) YAPP_Malloc((size_t) iNFFT,
                                          sizeof(float),
                                          YAPP_TRUE);
      if (NULL == g_pfSumPowX)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
      g_pfSumPowY = (float *) YAPP_Malloc((size_t) iNFFT,
                                          sizeof(float),
                                          YAPP_TRUE);
      if (NULL == g_pfSumPowY)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
      g_pfSumStokesRe = (float *) YAPP_Malloc((size_t) iNFFT,
                                              sizeof(float),
                                              YAPP_TRUE);
      if (NULL == g_pfSumStokesRe)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
      g_pfSumStokesIm = (float *) YAPP_Malloc((size_t) iNFFT,
                                              sizeof(float),
                                              YAPP_TRUE);
      if (NULL == g_pfSumStokesIm)
      {
          (void) fprintf(stderr,
                         "ERROR: Memory allocation failed! %s.\n",
                         strerror(errno));
          return YAPP_RET_ERROR;
      }
  
      /* create plans */
      g_stPlanX = fftwf_plan_dft_1d(iNFFT,
                                    g_pfcFFTInX,
                                    g_pfcFFTOutX,
                                    FFTW_FORWARD,
                                    FFTW_MEASURE);
      g_stPlanY = fftwf_plan_dft_1d(iNFFT,
                                    g_pfcFFTInY,
                                    g_pfcFFTOutY,
                                    FFTW_FORWARD,
                                    FFTW_MEASURE);
  
      return YAPP_RET_SUCCESS;
  }
  
  /* function that reads data from input buffer */
  int ReadData(char cIsFirst, int iNTaps, int iNFFT)
  {
      int iRet = YAPP_RET_SUCCESS;
      int i = 0;
      int j = 0;
      int k = 0;
  
      if (cIsFirst)
      {
          /* read first iNTaps blocks of data */
          /* ASSUMPTION: there is at least iNTaps blocks of data */
          for (i = 0; i < iNTaps; ++i)
          {
              iRet = read(g_iFileData, g_astPFBData[i].pcData, iNFFT * NUM_BYTES_PER_SAMP);
              if (iRet < YAPP_RET_SUCCESS)
              {
                  (void) fprintf(stderr,
                                 "ERROR: Data reading failed! %s.\n",
                                 strerror(errno));
                  return YAPP_RET_ERROR;
              }
  
              /* unpack data */
              /* assuming real and imaginary parts are interleaved, and X and Y are
                 interleaved, like so:
                 reX, imX, reY, imY, ... */
              j = 0;
              for (k = 0; k < NUM_BYTES_PER_SAMP * iNFFT; k += NUM_BYTES_PER_SAMP)
              {
                  g_astPFBData[i].pfcDataX[j][0]
                                          = (float) g_astPFBData[i].pcData[k];
                  g_astPFBData[i].pfcDataX[j][1]
                                          = (float) g_astPFBData[i].pcData[k+1];
                  g_astPFBData[i].pfcDataY[j][0]
                                          = (float) g_astPFBData[i].pcData[k+2];
                  g_astPFBData[i].pfcDataY[j][1]
                                          = (float) g_astPFBData[i].pcData[k+3];
                  ++j;
              }
          }
          g_iPFBWriteIdx = 0;     /* next write into the first buffer */
          g_iPFBReadIdx = 0;      /* PFB to be performed from first buffer */
      }
      else
      {
          /* write new data to the write buffer */
          iRet = read(g_iFileData, g_astPFBData[g_iPFBWriteIdx].pcData, iNFFT * NUM_BYTES_PER_SAMP);
          if (iRet < YAPP_RET_SUCCESS)
          {
              (void) fprintf(stderr,
                             "ERROR: Data reading failed! %s.\n",
                             strerror(errno));
              return YAPP_RET_ERROR;
          }
          if (iRet < iNFFT * NUM_BYTES_PER_SAMP)
          {
              g_iPFBReadIdx = g_astPFBData[g_iPFBReadIdx].iNextIdx;
              return YAPP_RET_READDONE;
          }
  
          /* unpack data */
          /* assuming real and imaginary parts are interleaved, and X and Y are
             interleaved, like so:
             Re(X), Im(X), Re(Y), Im(Y), ... */
          j = 0;
          for (i = 0; i < NUM_BYTES_PER_SAMP * iNFFT; i += NUM_BYTES_PER_SAMP)
          {
              g_astPFBData[g_iPFBWriteIdx].pfcDataX[j][0]
                                  = (float) g_astPFBData[g_iPFBWriteIdx].pcData[i];
              g_astPFBData[g_iPFBWriteIdx].pfcDataX[j][1]
                                  = (float) g_astPFBData[g_iPFBWriteIdx].pcData[i+1];
              g_astPFBData[g_iPFBWriteIdx].pfcDataY[j][0]
                                  = (float) g_astPFBData[g_iPFBWriteIdx].pcData[i+2];
              g_astPFBData[g_iPFBWriteIdx].pfcDataY[j][1]
                                  = (float) g_astPFBData[g_iPFBWriteIdx].pcData[i+3];
              ++j;
          }
          g_iPFBWriteIdx = g_astPFBData[g_iPFBWriteIdx].iNextIdx;
          g_iPFBReadIdx = g_astPFBData[g_iPFBReadIdx].iNextIdx;
      }
  
      return YAPP_RET_SUCCESS;
  }
  
  /* function that performs the PFB */
  int DoPFB(int iNTaps, int iNFFT)
  {
      int i = 0;
      int j = 0;
      int k = 0;
      int iCoeffStartIdx = 0;
  
      /* reset memory */
      (void) memset(g_pfcFFTInX, '\0', iNFFT * sizeof(fftwf_complex));
      (void) memset(g_pfcFFTInY, '\0', iNFFT * sizeof(fftwf_complex));
  
      i = g_iPFBReadIdx;
      for (j = 0; j < iNTaps; ++j)
      {
          iCoeffStartIdx = j * iNFFT;
          for (k = 0; k < iNFFT; ++k)
          {
              g_pfcFFTInX[k][0] += g_astPFBData[i].pfcDataX[k][0] * g_pfPFBCoeff[iCoeffStartIdx+k];
              g_pfcFFTInX[k][1] += g_astPFBData[i].pfcDataX[k][1] * g_pfPFBCoeff[iCoeffStartIdx+k];
              g_pfcFFTInY[k][0] += g_astPFBData[i].pfcDataY[k][0] * g_pfPFBCoeff[iCoeffStartIdx+k];
              g_pfcFFTInY[k][1] += g_astPFBData[i].pfcDataY[k][1] * g_pfPFBCoeff[iCoeffStartIdx+k];
          }
          i = g_astPFBData[i].iNextIdx;
      }
  
      return YAPP_RET_SUCCESS;
  }
  
  int CopyDataForFFT(int iNFFT)
  {
      (void) memcpy(g_pfcFFTInX,
                    g_astPFBData[g_iPFBReadIdx].pfcDataX,
                    iNFFT * sizeof(fftwf_complex));
      (void) memcpy(g_pfcFFTInY,
                    g_astPFBData[g_iPFBReadIdx].pfcDataY,
                    iNFFT * sizeof(fftwf_complex));
  
      return YAPP_RET_SUCCESS;
  }
  
  /* function that performs the FFT */
  int DoFFT()
  {
      /* execute plan */
      fftwf_execute(g_stPlanX);
      fftwf_execute(g_stPlanY);
  
      return YAPP_RET_SUCCESS;
  }
  
  /* function that frees resources */
  void CleanUp(int iNTaps)
  {
      int i = 0;
  
      /* free resources */
      for (i = 0; i < iNTaps; ++i)
      {
          if (g_astPFBData[i].pfcDataX != NULL)
          {
              fftwf_free(g_astPFBData[i].pfcDataX);
              g_astPFBData[i].pfcDataX = NULL;
          }
          if (g_astPFBData[i].pfcDataY != NULL)
          {
              fftwf_free(g_astPFBData[i].pfcDataY);
              g_astPFBData[i].pfcDataY = NULL;
          }
      }
      if (g_pfcFFTInX != NULL)
      {
          fftwf_free(g_pfcFFTInX);
          g_pfcFFTInX = NULL;
      }
      if (g_pfcFFTInY != NULL)
      {
          fftwf_free(g_pfcFFTInY);
          g_pfcFFTInY = NULL;
      }
      if (g_pfcFFTOutX != NULL)
      {
          fftwf_free(g_pfcFFTOutX);
          g_pfcFFTOutX = NULL;
      }
      if (g_pfcFFTOutY != NULL)
      {
          fftwf_free(g_pfcFFTOutY);
          g_pfcFFTOutY = NULL;
      }
  
      /* destroy plans */
      fftwf_destroy_plan(g_stPlanX);
      fftwf_destroy_plan(g_stPlanY);
  
      fftwf_cleanup();
  
      if (g_iFileData != 0)
      {
          (void) close(g_iFileData);
          g_iFileData = 0;
      }
  
      /* free all memory allocated using YAPP_Malloc() */
      YAPP_CleanUp();
  
      return;
  }
  
  /*
   * Prints usage information
   */
  void PrintUsage(const char *pcProgName)
  {
      (void) printf("Usage: %s [options] <data-file>\n",
                    pcProgName);
      (void) printf("    -h  --help                          ");
      (void) printf("Display this usage information\n");
      (void) printf("    -s  --skip <time>                   ");
      (void) printf("The length of data in seconds, to be\n");
      (void) printf("                                        ");
      (void) printf("skipped\n");
      (void) printf("                                        ");
      (void) printf("(default is 0 s)\n");
      (void) printf("    -p  --proc <time>                   ");
      (void) printf("The length of data in seconds, to be\n");
      (void) printf("                                        ");
      (void) printf("processed\n");
      (void) printf("                                        ");
      (void) printf("(default is all)\n");
      (void) printf("    -l  --tsamp <tsamp>                 ");
      (void) printf("Sampling time in s\n");
      (void) printf("    -f  --centre-freq <freq>            ");
      (void) printf("Centre frequency of observing band,\n");
      (void) printf("                                        ");
      (void) printf("in MHz\n");
      (void) printf("    -b  --pfb                           ");
      (void) printf("Do polyphase filter bank\n");
      (void) printf("    -t  --ntaps <ntaps>                 ");
      (void) printf("Number of taps in the PFB\n");
      (void) printf("                                        ");
      (void) printf("(default is 8)\n");
      (void) printf("    -n  --nfft <nfft>                   ");
      (void) printf("Number of points in FFT\n");
      (void) printf("                                        ");
      (void) printf("(default is 1024)\n");
      (void) printf("    -a  --nacc <acc-len>                ");
      (void) printf("Number of spectra to accumulate\n");
      (void) printf("                                        ");
      (void) printf("(default is 1)\n");
      (void) printf("    -o  --obs-site                      ");
      (void) printf("Observatory name\n");
      (void) printf("    -r  --pulsar                        ");
      (void) printf("Pulsar name\n");
      (void) printf("    -g  --graphics                      ");
      (void) printf("Turn on plotting\n");
      (void) printf("    -i  --invert                        ");
      (void) printf("Invert background and foreground\n");
      (void) printf("                                        ");
      (void) printf("colours in plots\n");
      (void) printf("    -e  --non-interactive               ");
      (void) printf("Run in non-interactive mode\n");
      (void) printf("    -v  --version                       ");
      (void) printf("Display the version\n");
  
      return;
  }
  
  