Changeset 2177


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
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/python/scantable.py

    r2150 r2177  
    10831083        """\
    10841084        Flag the data in 'lag' space by providing a frequency to remove.
    1085         Flagged data in the scantable gets interpolated over the region.
     1085        Flagged data in the scantable get interpolated over the region.
    10861086        No taper is applied.
    10871087
    … …  
    11181118        else:
    11191119            return s
     1120
     1121    @asaplog_post_dec
     1122    def fft(self, getrealimag=False, rowno=[]):
     1123        """\
     1124        Apply FFT to the spectra.
     1125        Flagged data in the scantable get interpolated over the region.
     1126
     1127        Parameters:
     1128       
     1129            getrealimag:    If True, returns the real and imaginary part
     1130                            values of the complex results.
     1131                            If False (the default), returns the amplitude
     1132                            (absolute value) normalised with Ndata/2 and
     1133                            phase (argument, in unit of radian).
     1134
     1135            rowno:          The row number(s) to be processed. int, list
     1136                            and tuple are accepted. By default ([]), FFT
     1137                            is applied to the whole data.
     1138
     1139        Returns:
     1140
     1141            A dictionary containing two keys, i.e., 'real' and 'imag' for
     1142            getrealimag = True, or 'ampl' and 'phase' for getrealimag = False,
     1143            respectively.
     1144            The value for each key is a list of lists containing the FFT
     1145            results from each spectrum.
     1146           
     1147        """
     1148        if getrealimag:
     1149            res = {"real":[], "imag":[]}  # return real and imaginary values
     1150        else:
     1151            res = {"ampl":[], "phase":[]} # return amplitude and phase(argument)
     1152       
     1153        if isinstance(rowno, int):
     1154            rowno = [rowno]
     1155        elif not (isinstance(rowno, list) or isinstance(rowno, tuple)):
     1156            raise TypeError("The row number(s) should be int, list or tuple.")
     1157       
     1158        nrow = len(rowno)
     1159        if nrow == 0: nrow = self.nrow()  # by default, calculate for all rows.
     1160       
     1161        fspec = self._math._fft(self, rowno, getrealimag)
     1162        nspec = len(fspec)/nrow
     1163       
     1164        i = 0
     1165        while (i < nrow):
     1166            v1 = []
     1167            v2 = []
     1168           
     1169            j = 0
     1170            while (j < nspec):
     1171                k = i*nspec + j
     1172                v1.append(fspec[k])
     1173                v2.append(fspec[k+1])
     1174                j += 2
     1175           
     1176            if getrealimag:
     1177                res["real"].append(v1)
     1178                res["imag"].append(v2)
     1179            else:
     1180                res["ampl"].append(v1)
     1181                res["phase"].append(v2)
     1182           
     1183            i += 1
     1184
     1185        return res
    11201186
    11211187    @asaplog_post_dec

  • 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.