hilbert.cpp

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 3.0 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
//  
//  Copyright 2004, The University of Texas at Austin
//
//============================================================================

//============================================================================
//
//This software developed by Applied Research Laboratories at the University of
//Texas at Austin, under contract to an agency or agencies within the U.S. 
//Department of Defense. The U.S. Government retains all rights to use,
//duplicate, distribute, disclose, or release this software. 
//
//Pursuant to DoD Directive 523024 
//
// DISTRIBUTION STATEMENT A: This software has been approved for public 
//                           release, distribution is unlimited.
//
//=============================================================================

/*

Program to perform Hilbert transform on samples from output file generated by SiGe SE4110L.  The SiGe chip outputs real 2-bit samples.  The GPSCreations GPS1A board ships with a program "ogusb-lit60.exe" that outputs the samples in char format to a file.  This program reads from such a file, performs the Hilbert transform, and outputs I and Q samples in an IQStream used by the GPStk.

The effective sampling rate is cut in half because N points of real data generate N/2 points of complex data.  The input sampling rate from the GPS1A is 16.368 MHz, so the output from this program will be 8.184 MHz.

Requires the FFTW library.

GCC: (not added to Jamfile yet, since this links to FFTW)
g++ -c -o hilbert.o -O -I. -I/.../gpstk/dev/apps/swrx -I/.../gpstk/dev/src hilbert.cpp

g++ -o hilbert hilbert.o /.../gpstk/dev/apps/swrx/simlib.a /.../gpstk/dev/src/libgpstk.a -lm -lstdc++ -lfftw3 -lm

*/

#include <iostream>
#include <fstream>
#include <stdio.h>
#include <math.h> 
#include <complex>
#include <fftw3.h>
#include <vector>
#include "CommandOption.hpp"
#include "CommandOptionParser.hpp"
#include "BasicFramework.hpp"
#include "IQStream.hpp"
using namespace gpstk;
using namespace std;
using namespace gpstk::StringUtils;

int main(int argc, char *argv[])
{
   CommandOptionWithAnyArg
      inputOpt('i', "input", "Location of input data file.",true);
   CommandOptionNoArg
      helpOption('h', "help", "Print usage. Repeat for more info. ");

   CommandOptionWithAnyArg
      timeOpt('t', "number-samples", "Number of C/A periods to input.");
   CommandOptionWithAnyArg
      quantizationOpt('q', "quantization",
                      "What type of IQ stream; 1, 2 or f. The default is 2.");
   
   string appDesc("Performs hilbert transform on GPS1A data.");
   CommandOptionParser cop(appDesc);
   cop.parseOptions(argc, argv);

   if (helpOption.getCount() || cop.hasErrors())
   {
      if (cop.hasErrors() && helpOption.getCount()==0)
      {
         cop.dumpErrors(cout);
         cout << "Use -h for help." << endl;
      }
      else
      {
         cop.displayUsage(cout);
      }
      exit(0);
   } 

   string f;
   f = inputOpt.getValue()[0].c_str();
   std::ifstream in;
   in.open(f.data(), std::ios::in | std::ios::binary);
   if(!in)
   {  cerr << "Could not open input file, exiting..." << endl;
      return 0;}

   IQStream *output;
   char quantization='2';   
   if (quantizationOpt.getCount())
      quantization = quantizationOpt.getValue()[0][0];
   switch (quantization)
   {
      case '1': output = new IQ1Stream(); break;
      case '2': output = new IQ2Stream(); break;
      case 'f': output = new IQFloatStream(); break;
      default:  output = new IQ2Stream(); break;
   }

   int numPeriods;
   if(timeOpt.getCount())
   {numPeriods=asInt(timeOpt.getValue().front());}
   // This is not currently implemented, right now we just read to
   // the end of the input file.

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
   int N = 64; // Default number of samples per transform.
   int N2 = N/2;
   char c[N];
   int count = 0;
   double *REAL;
   fftw_complex *IQ;
   fftw_complex *X;
   fftw_complex *X1;
   fftw_plan p,pb;
   
   using std::basic_ios;              // Set up IQStream
   output->copyfmt(std::cout);
   output->clear(std::cout.rdstate());
   output->basic_ios<char>::rdbuf(std::cout.rdbuf());
   output->filename = "<stdout>";

   REAL = (double*) malloc(sizeof(double) * N + 1); // input
   IQ = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N2); // final output
   X = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N); // intermediate
   X1 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N2);
   p = fftw_plan_dft_r2c_1d(N, REAL, X, FFTW_MEASURE); //FFT plan
   pb = fftw_plan_dft_1d(N2, X1, IQ, FFTW_BACKWARD, FFTW_MEASURE);

/*
   // Sloppy way to get samples from middle of file.
   char throwAway;
   for(int i = 0; i < 16368 * 2000; i++)
   {
      throwAway = fgetc(in);
   }
*/

   int size; // in bytes
   if(in.is_open())
   {
      in.seekg(0,ios::end);
      size = in.tellg();
      in.seekg(0, ios::beg);
   }
   while(size > 0)
   {
      // Read data into input array.
      in.read(c,N);
      size -= N;
      for(int i=0 ; i < N ; i++) 
      {
         REAL[i] = c[i];
      }

      // Execute FFT.
      fftw_execute(p);

      // Create new Hilbert components.
      for(int i=0; i < ( N/2 - 1 ) ; i++)
      {
         X1[i][0] = X[i][0];
         X1[i][1] = X[i][1];
      }
      
      // Execute Backward FFT.
      fftw_execute(pb); 

      for(int i=0 ; i < N2 ; i++)
      {  
         IQ[i][0] =   (  IQ[i][0] / N2 ) ;
         IQ[i][1] =   (  IQ[i][1] / N2 ) ;
      }

      // Output I and Q samples to IQStream.
      complex<float> o[N2];
      for(int i=0; i < (N2); i++)
      {
         o[i] = complex<float>(IQ[i][0],IQ[i][1]);
         *output << o[i];
      }
      count+=N;
   }//while

   fftw_destroy_plan(p);fftw_destroy_plan(pb);
   fftw_free(REAL); fftw_free(IQ); fftw_free(X); fftw_free(X1);
   return 0;
}