M3_format_otfs.c

/*******************************************************************************
*   M3:  M3_format_otfs.c                                                      *
*                                                                              *
*   Version 2.0 April 2007                                                     *
*                                                                              *
*   Copyright (C) 2004-2007  C.M. Cantalupo                                    *
*                                                                              *
*   This program is free software; you can redistribute it and/or modify       *
*   it under the terms of the GNU General Public License as published by       *
*   the Free Software Foundation; either version 2 of the License, or          *
*   (at your option) any later version.                                        *
*                                                                              *
*   This program 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 General Public License for more details.                               *
*                                                                              *
*   You should have received a copy of the GNU General Public License          *
*   along with this program; if not, write to the Free Software                *
*   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA  *
*                                                                              *
*******************************************************************************/

#include "M3.h"
#include "M3_internal.h"
#include "M3_format_otfs.h"
#include "M3_gcp.h"
#include <stdlib.h>
#include <stdio.h>

#include "rotateCoords.h"
#define M3_CONTENT_LENGTH 256


void M3_File_parseFormatString_otfs( M3_File *theFile, M3_FormatStruct_otfs *out );
void M3_File_readHeader_otfs( M3_File *theFile, void *data );
void M3_File_readData_otfs( M3_File *theFile, void *readArg, void *data );
void M3_File_readHeader_otfs_tod( M3_File *theFile, M3_TODheader *data );
void M3_File_readData_otfs_tod( M3_File *theFile, M3_Interval *readInterval, double *data );

void M3_CreateNoiseFilter( double fknee, double alpha, double fmin, double fsample, double NET, int64_t corr, double *inv_tt_noise );


void M3_File_parseFormatString_otfs( M3_File *theFile, M3_FormatStruct_otfs *out )
{
  M3_FormatStruct_otfs formatStruct;
  char *s1, *s2;
  char errorString[M3_CONTENT_LENGTH];
  static char lastFormatString[M3_CONTENT_LENGTH] = {0};
  static M3_FormatStruct_otfs lastFormatStruct;

   /* Build in a last call cache, as we will often be parsing the same format strings over and over again.  */
  if( strcmp( theFile->format, lastFormatString ) == 0 )
  {
    *out = lastFormatStruct;
    return;
  }
  
  /* Set default values */
  out->prngSeed = 0;
  out->prngFirstBlock = 0;
  out->prngNumBlock = 4096;
  out->prngOffset = 0;
  out->prngSigma = 1;
  out->firstSample = 0;
  out->lastSample = 0;
  out->oofSimFormat = NULL;
  out->fknee = 0;
  out->alpha = 0;
  out->fmin = 0;
  out->fsample = 0;
  out->NET = 0;
  out->corLength = 0;

  sprintf( errorString, "Error interpreting format string: \"%s\"", theFile->format );

  /* Advance to the first parameter */
  s1 = strchr( theFile->format, '_');

  /* Loop over the parameters in the format string */
  while( s1 )
  {
    /* Skip '_' character */
    s1++;

    /* Switch over all of the parameters */
    if( strncmp( s1, "prngSeed", 8 ) == 0 )
    {
      s1 += 8;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->prngSeed = strtoll(s1, NULL, 10);
    }
    else if( strncmp( s1, "prngFirstBlock", 14 ) == 0 )
    {
      s1 += 14;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      out->prngFirstBlock = strtoll( s1, NULL, 10 );
    }
    else if( strncmp( s1, "prngNumBlock", 12 ) == 0 )
    {
      s1 += 12;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->prngNumBlock = strtoll( s1, NULL, 10 );
    }
    else if( strncmp( s1, "prngOffset", 10 ) == 0 )
    {
      s1 += 10;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->prngOffset = strtoll( s1, NULL, 10 );
    }
    else if( strncmp( s1, "prngSigma", 9 ) == 0 )
    {
      s1 += 9;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      out->prngSigma = atof( s1 );
    }
    else if( strncmp( s1, "firstSample", 11 ) == 0 )
    {
      s1 += 11;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->firstSample = strtoll(s1, NULL, 10 );
    }
    else if( strncmp( s1, "lastSample", 10 ) == 0 )
    {
      s1 += 10;
      M3_ErrorCheck( -1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->lastSample = strtoll(s1, NULL, 10 );
    }
    else if( strncmp( s1, "oofSimFormat", 12 ) == 0 )
    {
      s1 += 12;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT );
      s1++;
      out->oofSimFormat = s1;
    }
    else if( strncmp( s1, "fknee", 5 ) == 0 )
    {
      s1 += 5;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->fknee = atof(s1);
    }
    else if( strncmp( s1, "alpha", 5 ) == 0 )
    {
      s1 += 5;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->alpha = atof(s1);
    }
    else if( strncmp( s1, "fmin", 4 ) == 0 )
    {
      s1 += 4;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->fmin = atof(s1);
    }
    else if( strncmp( s1, "fsample", 7 ) == 0 )
    {
      s1 += 7;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->fsample = atof(s1);
    }
    else if( strncmp( s1, "NET", 3 ) == 0 )
    {
      s1 += 3;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->NET = atof(s1);
    }
    else if( strncmp( s1, "corLength", 9 ) == 0 )
    {
      s1 += 9;
      M3_ErrorCheck(-1, errorString, *s1 == '=', M3_EFLAG_DEFAULT);
      s1++;
      out->corLength = strtoll(s1, NULL, 10);
    }


    s1 = strchr( s1, '_');
  }


  /* Set last call cache */
  strcpy( lastFormatString, theFile->format );
  lastFormatStruct = *out;

}

void M3_File_readHeader_otfs_noise( M3_File *theFile, M3_NoiseHeader *header )
{
  M3_FormatStruct_otfs formatStruct;

  M3_File_parseFormatString_otfs( theFile, &formatStruct );

  header->firstSample = formatStruct.firstSample;
  header->lastSample = formatStruct.lastSample;
  header->corLength = formatStruct.corLength;

}


void M3_File_readData_otfs_noise( M3_File *theFile, double *noise )
{
  M3_FormatStruct_otfs formatStruct;

  M3_File_parseFormatString_otfs( theFile, &formatStruct );

  M3_CreateNoiseFilter( formatStruct.fknee, formatStruct.alpha, formatStruct.fmin, formatStruct.fsample, formatStruct.NET, formatStruct.corLength, noise ); 

}

void M3_File_readHeader_otfs( M3_File *theFile, void *data )
{
  /* Simple switch over sub-functions for each file type */
  switch( theFile->fileType )
  {
    case M3_TOD_FILE_TYPE:
      M3_File_readHeader_otfs_tod( theFile, (M3_TODheader *)data );
      break;
    case M3_NOISE_FILE_TYPE:
      M3_File_readHeader_otfs_noise(theFile, (M3_NoiseHeader *)data );
      break;  
    default:
      M3_ErrorCheck(-1, "M3 can only read otfs files for the tod and noise file type", 0, M3_EFLAG_DEFAULT );
      break;
  }
}

void M3_File_readData_otfs( M3_File *theFile, void *readArg, void *data )
{
  /* Simple switch over sub-functions for each file type */
  switch( theFile->fileType )
  {
    case M3_TOD_FILE_TYPE:
      M3_File_readData_otfs_tod( theFile, (M3_Interval*)readArg, (double *)data );
      break;
    case M3_NOISE_FILE_TYPE:
      M3_File_readData_otfs_noise( theFile, (double *)data );
      break;
    default:
      M3_ErrorCheck(-1, "M3 can only read otfs files for the tod noise file type", 0, M3_EFLAG_DEFAULT );
      break;
  }
}

void M3_File_readHeader_otfs_tod( M3_File *theFile, M3_TODheader *data )
{
  M3_FormatStruct_otfs formatStruct;
  char *strPtr;
  M3_NoiseHeader noiseHeader;

  M3_File_parseFormatString_otfs( theFile, &formatStruct);
  if( formatStruct.oofSimFormat )
  {
    theFile->fileType = M3_NOISE_FILE_TYPE;
    strPtr = theFile->format;
    theFile->format = formatStruct.oofSimFormat;
    M3_File_ReadHeader( theFile, &noiseHeader);
    data->firstSample = noiseHeader.firstSample;
    data->lastSample = noiseHeader.lastSample;
    theFile->param.noise.corLength = noiseHeader.corLength;
    theFile->format = strPtr;
    theFile->fileType = M3_TOD_FILE_TYPE;
  }
  else
  {
    data->firstSample = formatStruct.firstSample;
    data->lastSample = formatStruct.lastSample;
  }

}

void M3_File_readData_otfs_tod( M3_File *theFile, M3_Interval *readInterval, double *data )
{
  M3_FormatStruct_otfs formatStruct;
  int64_t offset, numReadEl, numFileEl, i, blockSize, j, blockIndex, blockOffset, thisNumEl, fftLength;
  double *noise, *filter, *forwardPlan, *backwardPlan, *nptr;
  double temp;
  int err;
  char *strPtr;
  M3_NoiseHeader noiseHeader;

  M3_File_parseFormatString_otfs( theFile, &formatStruct );

  numReadEl = readInterval->lastSample - readInterval->firstSample + 1;
  numFileEl = theFile->param.tod.interval.lastSample - theFile->param.tod.interval.firstSample + 1;

  if( !formatStruct.oofSimFormat )
  {
    offset = readInterval->firstSample - theFile->param.tod.interval.firstSample + formatStruct.prngOffset;
    
    blockSize = (numFileEl + formatStruct.prngOffset)/formatStruct.prngNumBlock + ((numFileEl + formatStruct.prngOffset)%formatStruct.prngNumBlock ? 1:0);
    blockIndex = offset/blockSize;
    blockOffset = offset%blockSize;

    thisNumEl = blockSize - blockOffset; 
    if( thisNumEl > numReadEl )
      thisNumEl = numReadEl;

    i = 0;
    /* Get data from part of the first used block */
    SRNG_Generate_WhiteNoise( formatStruct.prngSeed, formatStruct.prngFirstBlock + blockIndex, blockOffset, thisNumEl, formatStruct.prngSigma, data );
    i += thisNumEl;
    blockIndex++;
    
    /* Get data from the full blocks */
    while( numReadEl - i  > blockSize )
    {
      SRNG_Generate_WhiteNoise( formatStruct.prngSeed, formatStruct.prngFirstBlock + blockIndex, 0, blockSize, formatStruct.prngSigma, data + i );
      i += blockSize;
      blockIndex++;
    }
  
    /* Get the last part */
    if ( numReadEl - i )
    {
      thisNumEl = numReadEl - i;
      SRNG_Generate_WhiteNoise( formatStruct.prngSeed, formatStruct.prngFirstBlock + blockIndex, 0, thisNumEl, formatStruct.prngSigma, data + i );
    }
  }
  else
  {
#ifndef USE_LIBACML
    M3_ErrorCheck(-1, "Generation of noise requires ACML, FFTW not yet implemented", 0, M3_EFLAG_DEFAULT );
#else

    /* allocate and initizlize two arrays that are twice the length of the file */
    fftLength = pow(2,floor(log2(2*numFileEl-1))+1);

    filter = (double *)calloc(fftLength,sizeof(double));
    noise = (double *)calloc(fftLength,sizeof(double));
    forwardPlan = (double *)malloc((3*fftLength+100)*sizeof(double));
    backwardPlan = (double *)malloc((3*fftLength+100)*sizeof(double));
    M3_ErrorCheck(-1, "M3_File_readData_otfs_tod(): FFT vectors", (filter && noise && forwardPlan && backwardPlan ), M3_EFLAG_MALLOC_FSTRING );
    /* Initialize workspace */
    dzfft( 0, fftLength, noise, forwardPlan, &err );
    zdfft( 0, fftLength, noise, backwardPlan, &err );

    /* read in the inv_tt_noise file into the first bit of the filter vector.*/
    theFile->fileType = M3_NOISE_FILE_TYPE;
    strPtr = theFile->format;
    theFile->format = formatStruct.oofSimFormat;
    M3_File_ReadHeader( theFile, &noiseHeader);
    M3_ErrorCheck(-1, "Noise generation OTFS TOD has a correlation length that is longer than half the file length.", 2*noiseHeader.corLength < numFileEl, M3_EFLAG_DEFAULT);
    M3_File_ReadData( theFile, NULL, noise );
    theFile->format = strPtr;
    theFile->fileType = M3_TOD_FILE_TYPE;


    /* copy it into the last bit of the filter vector. */
    for( i = 1; i <= noiseHeader.corLength; i++)
      noise[fftLength - i] = noise[i];

    /* fft the filter vector */
    dzfft(1, fftLength, noise, forwardPlan, &err );

    /***********************************************************************
    *   Now we know that the Fourier domain filter is real valued and      *
    *   even, so convert the filter from Hermitian form to real (and       *
    *   even).  We also need to conjugate before we do the backward FFT    *
    *   of the filtered data, so change the sign of the imaginary part of  *
    *   the filter.                                                        *
    ***********************************************************************/

    /* invert the square root of each element to go from N^-1 -> N^1/2 
       (note that it is positive and real valued) */
    temp = pow(fftLength,-0.25);
    filter[0] = temp/sqrt(noise[0]);
    i = fftLength/2 + 1;
    for( j = 1; j < i; j++ )
    {
      filter[j] = temp/sqrt(noise[j]);
      filter[fftLength - j] = -1*filter[j];
    }

    /* fill the first half of the noise array with normal zero one noise */
    memset(noise, 0, sizeof(double)*fftLength);
    offset = formatStruct.prngOffset;
    numFileEl = theFile->param.tod.interval.lastSample - theFile->param.tod.interval.firstSample + 1;

    blockSize = (numFileEl + formatStruct.prngOffset)/formatStruct.prngNumBlock + ((numFileEl + formatStruct.prngOffset)%formatStruct.prngNumBlock ? 1:0);
    blockIndex = offset/blockSize;
    blockOffset = offset%blockSize;

    thisNumEl = blockSize - blockOffset; 
    if( thisNumEl > numFileEl )
      thisNumEl = numFileEl;

    i = 0;
    /* Get data from part of the first used block */
    SRNG_Generate_WhiteNoise( formatStruct.prngSeed, formatStruct.prngFirstBlock + blockIndex, blockOffset, thisNumEl, 1.0, noise );
    i += thisNumEl;
    blockIndex++;
    
    /* Get data from the full blocks */
    while( numFileEl - i  > blockSize )
    {
      SRNG_Generate_WhiteNoise( formatStruct.prngSeed, formatStruct.prngFirstBlock + blockIndex, 0, blockSize, 1.0, noise + i );
      i += blockSize;
      blockIndex++;
    }
  
    /* Get the last part */
    if ( numFileEl - i )
    {
      thisNumEl = numFileEl - i;
      SRNG_Generate_WhiteNoise( formatStruct.prngSeed, formatStruct.prngFirstBlock + blockIndex, 0, thisNumEl, 1.0, noise + i );
    }

    /* fft the noise vector */
    dzfft(1, fftLength, noise, forwardPlan, &err );     
    /* multiply the noise vector by the filter vector. */
    for( i = 0; i < fftLength; i++ )
      noise[i] *= filter[i];
    /* ifft the noise vector */
    zdfft(1, fftLength, noise, backwardPlan, &err );
    /* fill in the elements of the output array from the noise vector */
    nptr = noise + readInterval->firstSample - theFile->param.tod.interval.firstSample;
    for( i = 0; i < numReadEl; i++ )
      data[i] = nptr[i];
    free(backwardPlan);
    free(forwardPlan);
    free(noise);
    free(filter);
#endif

  }

}

---