Changeset 2177 for trunk/src


Ignore:
Timestamp:
05/20/11 18:53:58 (14 years ago)
Author:
WataruKawasaki
Message:

New Development: Yes

JIRA Issue: Yes CAS-2828

Ready for Test: Yes

Interface Changes:

What Interface Changed:

Test Programs:

Put in Release Notes:

Module(s): SD

Description: created a tool function sd.scantable.fft() to apply FFT for scantable data.


Location:
trunk/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/STMath.cpp

    r2163 r2177  
    1111//
    1212
     13#include <sstream>
     14
    1315#include <casa/iomanip.h>
    14 #include <casa/Exceptions/Error.h>
    15 #include <casa/Containers/Block.h>
    16 #include <casa/BasicSL/String.h>
    1716#include <casa/Arrays/MaskArrLogi.h>
    1817#include <casa/Arrays/MaskArrMath.h>
    … …  
    2120#include <casa/Arrays/Slice.h>
    2221#include <casa/Arrays/Slicer.h>
     22#include <casa/BasicSL/String.h>
     23#include <casa/Containers/Block.h>
    2324#include <casa/Containers/RecordField.h>
     25#include <casa/Exceptions/Error.h>
     26#include <casa/Logging/LogIO.h>
     27
     28#include <coordinates/Coordinates/CoordinateSystem.h>
     29#include <coordinates/Coordinates/CoordinateUtil.h>
     30#include <coordinates/Coordinates/FrequencyAligner.h>
     31#include <coordinates/Coordinates/SpectralCoordinate.h>
     32
     33#include <lattices/Lattices/LatticeUtilities.h>
     34
     35#include <scimath/Functionals/Polynomial.h>
     36#include <scimath/Mathematics/Convolver.h>
     37#include <scimath/Mathematics/VectorKernel.h>
     38
     39#include <tables/Tables/ExprNode.h>
     40#include <tables/Tables/ReadAsciiTable.h>
     41#include <tables/Tables/TableCopy.h>
     42#include <tables/Tables/TableIter.h>
     43#include <tables/Tables/TableParse.h>
     44#include <tables/Tables/TableRecord.h>
    2445#include <tables/Tables/TableRow.h>
    2546#include <tables/Tables/TableVector.h>
    2647#include <tables/Tables/TabVecMath.h>
    27 #include <tables/Tables/ExprNode.h>
    28 #include <tables/Tables/TableRecord.h>
    29 #include <tables/Tables/TableParse.h>
    30 #include <tables/Tables/ReadAsciiTable.h>
    31 #include <tables/Tables/TableIter.h>
    32 #include <tables/Tables/TableCopy.h>
    33 #include <scimath/Mathematics/FFTServer.h>
    34 
    35 #include <lattices/Lattices/LatticeUtilities.h>
    36 
    37 #include <coordinates/Coordinates/SpectralCoordinate.h>
    38 #include <coordinates/Coordinates/CoordinateSystem.h>
    39 #include <coordinates/Coordinates/CoordinateUtil.h>
    40 #include <coordinates/Coordinates/FrequencyAligner.h>
    41 
    42 #include <scimath/Mathematics/VectorKernel.h>
    43 #include <scimath/Mathematics/Convolver.h>
    44 #include <scimath/Functionals/Polynomial.h>
    4548
    4649#include <atnf/PKSIO/SrcType.h>
    47 
    48 #include <casa/Logging/LogIO.h>
    49 #include <sstream>
    5050
    5151#include "MathUtils.h"
    5252#include "RowAccumulator.h"
    5353#include "STAttr.h"
     54#include "STMath.h"
    5455#include "STSelector.h"
    5556
    56 #include "STMath.h"
    5757using namespace casa;
    58 
    5958using namespace asap;
    6059
    … …  
    28212820}
    28222821
    2823 CountedPtr< Scantable >
    2824   asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
     2822std::vector<float>
     2823  asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
     2824                     const std::vector<int>& whichrow,
     2825                     bool getRealImag )
     2826{
     2827  std::vector<float> res;
     2828  Table tab = in->table();
     2829  ArrayColumn<Float> specCol(tab, "SPECTRA");
     2830  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
     2831  FFTServer<Float,Complex> ffts;
     2832
     2833  if (whichrow.size() < 1) {          // for all rows (by default)
     2834    int nrow = int(tab.nrow());
     2835    for (int i = 0; i < nrow; ++i) {
     2836      Vector<Float> spec = specCol(i);
     2837      Vector<uChar> flag = flagCol(i);
     2838      doFFT(res, ffts, getRealImag, spec, flag);
     2839    }
     2840  } else {                           // for specified rows
     2841    for (uInt i = 0; i < whichrow.size(); ++i) {
     2842      Vector<Float> spec = specCol(whichrow[i]);
     2843      Vector<uChar> flag = flagCol(whichrow[i]);
     2844      doFFT(res, ffts, getRealImag, spec, flag);
     2845    }
     2846  }
     2847
     2848  return res;
     2849}
     2850
     2851void asap::STMath::doFFT( std::vector<float>& out,
     2852                          FFTServer<Float,Complex>& ffts,
     2853                          bool getRealImag,
     2854                          Vector<Float>& spec,
     2855                          Vector<uChar>& flag )
     2856{
     2857  doZeroOrderInterpolation(spec, flag);
     2858
     2859  Vector<Complex> fftout;
     2860  ffts.fft0(fftout, spec);
     2861
     2862  float norm = float(2.0/double(spec.size()));
     2863  if (getRealImag) {
     2864    for (uInt j = 0; j < fftout.size(); ++j) {
     2865      out.push_back(real(fftout[j])*norm);
     2866      out.push_back(imag(fftout[j])*norm);
     2867    }
     2868  } else {
     2869    for (uInt j = 0; j < fftout.size(); ++j) {
     2870      out.push_back(abs(fftout[j])*norm);
     2871      out.push_back(arg(fftout[j]));
     2872    }
     2873  }
     2874
     2875}
     2876
     2877void asap::STMath::doZeroOrderInterpolation( Vector<Float>& spec,
     2878                                             Vector<uChar>& flag )
     2879{
     2880  int fstart = -1;
     2881  int fend = -1;
     2882  for (uInt i = 0; i < flag.nelements(); ++i) {
     2883    if (flag[i] > 0) {
     2884      fstart = i;
     2885      while (flag[i] > 0 && i < flag.nelements()) {
     2886        fend = i;
     2887        i++;
     2888      }
     2889    }
     2890
     2891    // execute interpolation as the following criteria:
     2892    // (1) for a masked region inside the spectrum, replace the spectral
     2893    //     values with the mean of those at the two channels just outside
     2894    //     the both edges of the masked region.
     2895    // (2) for a masked region at the spectral edge, replace the values
     2896    //     with the one at the nearest non-masked channel.
     2897    //     (ZOH, but bilateral)
     2898    Float interp = 0.0;
     2899    if (fstart-1 > 0) {
     2900      interp = spec[fstart-1];
     2901      if (fend+1 < Int(spec.nelements())) {
     2902        interp = (interp + spec[fend+1]) / 2.0;
     2903      }
     2904    } else {
     2905      interp = spec[fend+1];
     2906    }
     2907    if (fstart > -1 && fend > -1) {
     2908      for (int j = fstart; j <= fend; ++j) {
     2909        spec[j] = interp;
     2910      }
     2911    }
     2912
     2913    fstart = -1;
     2914    fend = -1;
     2915  }
     2916}
     2917
     2918CountedPtr<Scantable>
     2919  asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
    28252920                         double start, double end,
    2826                          const std::string& mode)
    2827 {
    2828   CountedPtr< Scantable > out = getScantable(in, false);
     2921                         const std::string& mode )
     2922{
     2923  CountedPtr<Scantable> out = getScantable(in, false);
    28292924  Table& tout = out->table();
    28302925  TableIterator iter(tout, "FREQ_ID");
    28312926  FFTServer<Float,Complex> ffts;
     2927
    28322928  while ( !iter.pastEnd() ) {
    28332929    Table tab = iter.table();
    … …  
    28392935    ArrayColumn<Float> specCol(tab, "SPECTRA");
    28402936    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
     2937
    28412938    for (int i=0; i<int(tab.nrow()); ++i) {
    28422939      Vector<Float> spec = specCol(i);
    28432940      Vector<uChar> flag = flagCol(i);
    2844       int fstart = -1;
    2845       int fend = -1;
    2846       for (unsigned int k=0; k < flag.nelements(); ++k ) {
    2847         if (flag[k] > 0) {
    2848           fstart = k;
    2849           while (flag[k] > 0 && k < flag.nelements()) {
    2850             fend = k;
    2851             k++;
    2852           }
    2853         }
    2854         Float interp = 0.0;
    2855         if (fstart-1 > 0 ) {
    2856           interp = spec[fstart-1];
    2857           if (fend+1 < Int(spec.nelements())) {
    2858             interp = (interp+spec[fend+1])/2.0;
    2859           }
    2860         } else {
    2861           interp = spec[fend+1];
    2862         }
    2863         if (fstart > -1 && fend > -1) {
    2864           for (int j=fstart;j<=fend;++j) {
    2865             spec[j] = interp;
    2866           }
    2867         }
    2868         fstart =-1;
    2869         fend = -1;
    2870       }
     2941      doZeroOrderInterpolation(spec, flag);
     2942
    28712943      Vector<Complex> lags;
    28722944      ffts.fft0(lags, spec);
     2945
    28732946      Int lag0(start+0.5);
    28742947      Int lag1(end+0.5);
    … …  
    28782951      }
    28792952      Int lstart =  max(0, lag0);
    2880       Int lend =  min(Int(lags.nelements()-1), lag1);
     2953      Int lend   =  min(Int(lags.nelements()-1), lag1);
    28812954      if (lstart == lend) {
    28822955        lags[lstart] = Complex(0.0);
    … …  
    28912964        }
    28922965      }
     2966
    28932967      ffts.fft0(spec, lags);
     2968
    28942969      specCol.put(i, spec);
    28952970    }

  • trunk/src/STMath.h

    r1937 r2177  
    1313#define ASAPSTMATH_H
    1414
     15#include <map>
    1516#include <string>
    16 #include <map>
    1717
    1818#include <casa/aips.h>
     19#include <casa/Arrays/Vector.h>
     20#include <casa/BasicSL/String.h>
    1921#include <casa/Utilities/CountedPtr.h>
    20 #include <casa/BasicSL/String.h>
    21 #include <casa/Arrays/Vector.h>
     22
     23#include <scimath/Mathematics/FFTServer.h>
    2224#include <scimath/Mathematics/InterpolateArray1D.h>
    2325
     26#include "Logger.h"
    2427#include "Scantable.h"
    2528#include "STDefs.h"
    2629#include "STPol.h"
    27 #include "Logger.h"
    2830
    2931namespace asap {
    … …  
    345347    lagFlag( const casa::CountedPtr<Scantable>& in, double start,
    346348             double end, const std::string& mode="frequency");
     349
     350  std::vector<float>
     351    fft( const casa::CountedPtr<Scantable>& in,
     352         const std::vector<int>& whichrow,
     353         bool getRealImag=false );
    347354
    348355  // test for average spectra with different channel/resolution
    … …  
    426433                                        int index ) ;
    427434  double getMJD( string strtime ) ;
     435  void doFFT( std::vector<float>& out,
     436              casa::FFTServer< casa::Float, casa::Complex >& ffts,
     437              bool getRealImag,
     438              casa::Vector<casa::Float>& spec,
     439              casa::Vector<casa::uChar>& flag ) ;
     440  void doZeroOrderInterpolation( casa::Vector<casa::Float>& spec,
     441                                 casa::Vector<casa::uChar>& flag) ;
    428442
    429443  bool insitu_;

  • trunk/src/STMathWrapper.h

    r1907 r2177  
    213213                                            mode)); }
    214214
     215  std::vector<float> fft( const ScantableWrapper& in,
     216                          const std::vector<int>& whichrow,
     217                          const bool getRealImag=false )
     218  { return STMath::fft(in.getCP(), whichrow, getRealImag); }
     219
    215220  // test for average spectra with different channel/resolution
    216221  ScantableWrapper

  • trunk/src/python_STMath.cpp

    r1907 r2177  
    7777        .def("_mx_extract", &STMathWrapper::mxExtract)
    7878        .def("_lag_flag", &STMathWrapper::lagFlag)
     79        .def("_fft", &STMathWrapper::fft)
    7980        // testing average spectra with different channel/resolution
    8081        .def("_new_average", &STMathWrapper::new_average)
Note: See TracChangeset for help on using the changeset viewer.