Changeset 2177

Timestamp:

05/20/11 18:53:58 (13 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:

• python/scantable.py (modified) (2 diffs)
• src/STMath.cpp (modified) (6 diffs)
• src/STMath.h (modified) (3 diffs)
• src/STMathWrapper.h (modified) (1 diff)
• src/python_STMath.cpp (modified) (1 diff)

trunk/python/scantable.py

@@ -1083 +1083 @@
         """\
         Flag the data in 'lag' space by providing a frequency to remove.
-        Flagged data in the scantable gets interpolated over the region.
+        Flagged data in the scantable get interpolated over the region.
         No taper is applied.
 
@@ -1118 +1118 @@
         else:
             return s
+
+    @asaplog_post_dec
+    def fft(self, getrealimag=False, rowno=[]):
+        """\
+        Apply FFT to the spectra.
+        Flagged data in the scantable get interpolated over the region.
+
+        Parameters:
+        
+            getrealimag:    If True, returns the real and imaginary part 
+                            values of the complex results.
+                            If False (the default), returns the amplitude
+                            (absolute value) normalised with Ndata/2 and
+                            phase (argument, in unit of radian).
+
+            rowno:          The row number(s) to be processed. int, list 
+                            and tuple are accepted. By default ([]), FFT 
+                            is applied to the whole data.
+
+        Returns:
+
+            A dictionary containing two keys, i.e., 'real' and 'imag' for
+            getrealimag = True, or 'ampl' and 'phase' for getrealimag = False,
+            respectively.
+            The value for each key is a list of lists containing the FFT
+            results from each spectrum.
+            
+        """
+        if getrealimag:
+            res = {"real":[], "imag":[]}  # return real and imaginary values
+        else:
+            res = {"ampl":[], "phase":[]} # return amplitude and phase(argument)
+        
+        if isinstance(rowno, int):
+            rowno = [rowno]
+        elif not (isinstance(rowno, list) or isinstance(rowno, tuple)):
+            raise TypeError("The row number(s) should be int, list or tuple.")
+        
+        nrow = len(rowno)
+        if nrow == 0: nrow = self.nrow()  # by default, calculate for all rows. 
+        
+        fspec = self._math._fft(self, rowno, getrealimag)
+        nspec = len(fspec)/nrow
+        
+        i = 0
+        while (i < nrow):
+            v1 = []
+            v2 = []
+            
+            j = 0
+            while (j < nspec):
+                k = i*nspec + j
+                v1.append(fspec[k])
+                v2.append(fspec[k+1])
+                j += 2
+            
+            if getrealimag:
+                res["real"].append(v1)
+                res["imag"].append(v2)
+            else:
+                res["ampl"].append(v1)
+                res["phase"].append(v2)
+            
+            i += 1
+
+        return res
 
     @asaplog_post_dec

trunk/src/STMath.cpp

@@ -11 +11 @@
 //
 
+#include <sstream>
+
 #include <casa/iomanip.h>
-#include <casa/Exceptions/Error.h>
-#include <casa/Containers/Block.h>
-#include <casa/BasicSL/String.h>
 #include <casa/Arrays/MaskArrLogi.h>
 #include <casa/Arrays/MaskArrMath.h>
@@ -21 +20 @@
 #include <casa/Arrays/Slice.h>
 #include <casa/Arrays/Slicer.h>
+#include <casa/BasicSL/String.h>
+#include <casa/Containers/Block.h>
 #include <casa/Containers/RecordField.h>
+#include <casa/Exceptions/Error.h>
+#include <casa/Logging/LogIO.h>
+
+#include <coordinates/Coordinates/CoordinateSystem.h>
+#include <coordinates/Coordinates/CoordinateUtil.h>
+#include <coordinates/Coordinates/FrequencyAligner.h>
+#include <coordinates/Coordinates/SpectralCoordinate.h>
+
+#include <lattices/Lattices/LatticeUtilities.h>
+
+#include <scimath/Functionals/Polynomial.h>
+#include <scimath/Mathematics/Convolver.h>
+#include <scimath/Mathematics/VectorKernel.h>
+
+#include <tables/Tables/ExprNode.h>
+#include <tables/Tables/ReadAsciiTable.h>
+#include <tables/Tables/TableCopy.h>
+#include <tables/Tables/TableIter.h>
+#include <tables/Tables/TableParse.h>
+#include <tables/Tables/TableRecord.h>
 #include <tables/Tables/TableRow.h>
 #include <tables/Tables/TableVector.h>
 #include <tables/Tables/TabVecMath.h>
-#include <tables/Tables/ExprNode.h>
-#include <tables/Tables/TableRecord.h>
-#include <tables/Tables/TableParse.h>
-#include <tables/Tables/ReadAsciiTable.h>
-#include <tables/Tables/TableIter.h>
-#include <tables/Tables/TableCopy.h>
-#include <scimath/Mathematics/FFTServer.h>
-
-#include <lattices/Lattices/LatticeUtilities.h>
-
-#include <coordinates/Coordinates/SpectralCoordinate.h>
-#include <coordinates/Coordinates/CoordinateSystem.h>
-#include <coordinates/Coordinates/CoordinateUtil.h>
-#include <coordinates/Coordinates/FrequencyAligner.h>
-
-#include <scimath/Mathematics/VectorKernel.h>
-#include <scimath/Mathematics/Convolver.h>
-#include <scimath/Functionals/Polynomial.h>
 
 #include <atnf/PKSIO/SrcType.h>
-
-#include <casa/Logging/LogIO.h>
-#include <sstream>
 
 #include "MathUtils.h"
 #include "RowAccumulator.h"
 #include "STAttr.h"
+#include "STMath.h"
 #include "STSelector.h"
 
-#include "STMath.h"
 using namespace casa;
-
 using namespace asap;
 
@@ -2821 +2820 @@
 }
 
-CountedPtr< Scantable >
-  asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
+std::vector<float>
+  asap::STMath::fft( const casa::CountedPtr< Scantable > & in, 
+                     const std::vector<int>& whichrow, 
+                     bool getRealImag )
+{
+  std::vector<float> res;
+  Table tab = in->table();
+  ArrayColumn<Float> specCol(tab, "SPECTRA");
+  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+  FFTServer<Float,Complex> ffts;
+
+  if (whichrow.size() < 1) {          // for all rows (by default)
+    int nrow = int(tab.nrow());
+    for (int i = 0; i < nrow; ++i) {
+      Vector<Float> spec = specCol(i);
+      Vector<uChar> flag = flagCol(i);
+      doFFT(res, ffts, getRealImag, spec, flag);
+    }
+  } else {                           // for specified rows
+    for (uInt i = 0; i < whichrow.size(); ++i) {
+      Vector<Float> spec = specCol(whichrow[i]);
+      Vector<uChar> flag = flagCol(whichrow[i]);
+      doFFT(res, ffts, getRealImag, spec, flag);
+    }
+  }
+
+  return res;
+}
+
+void asap::STMath::doFFT( std::vector<float>& out, 
+                          FFTServer<Float,Complex>& ffts, 
+                          bool getRealImag,
+                          Vector<Float>& spec, 
+                          Vector<uChar>& flag )
+{
+  doZeroOrderInterpolation(spec, flag);
+
+  Vector<Complex> fftout;
+  ffts.fft0(fftout, spec);
+
+  float norm = float(2.0/double(spec.size()));
+  if (getRealImag) {
+    for (uInt j = 0; j < fftout.size(); ++j) {
+      out.push_back(real(fftout[j])*norm);
+      out.push_back(imag(fftout[j])*norm);
+    }
+  } else {
+    for (uInt j = 0; j < fftout.size(); ++j) {
+      out.push_back(abs(fftout[j])*norm);
+      out.push_back(arg(fftout[j]));
+    }
+  }
+
+}
+
+void asap::STMath::doZeroOrderInterpolation( Vector<Float>& spec, 
+                                             Vector<uChar>& flag )
+{
+  int fstart = -1;
+  int fend = -1;
+  for (uInt i = 0; i < flag.nelements(); ++i) {
+    if (flag[i] > 0) {
+      fstart = i;
+      while (flag[i] > 0 && i < flag.nelements()) {
+        fend = i;
+        i++;
+      }
+    }
+
+    // execute interpolation as the following criteria:
+    // (1) for a masked region inside the spectrum, replace the spectral 
+    //     values with the mean of those at the two channels just outside 
+    //     the both edges of the masked region. 
+    // (2) for a masked region at the spectral edge, replace the values 
+    //     with the one at the nearest non-masked channel. 
+    //     (ZOH, but bilateral)
+    Float interp = 0.0;
+    if (fstart-1 > 0) {
+      interp = spec[fstart-1];
+      if (fend+1 < Int(spec.nelements())) {
+        interp = (interp + spec[fend+1]) / 2.0;
+      }
+    } else {
+      interp = spec[fend+1];
+    }
+    if (fstart > -1 && fend > -1) {
+      for (int j = fstart; j <= fend; ++j) {
+        spec[j] = interp;
+      }
+    }
+
+    fstart = -1;
+    fend = -1;
+  }
+}
+
+CountedPtr<Scantable>
+  asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
                          double start, double end,
-                         const std::string& mode)
-{
-  CountedPtr< Scantable > out = getScantable(in, false);
+                         const std::string& mode )
+{
+  CountedPtr<Scantable> out = getScantable(in, false);
   Table& tout = out->table();
   TableIterator iter(tout, "FREQ_ID");
   FFTServer<Float,Complex> ffts;
+
   while ( !iter.pastEnd() ) {
     Table tab = iter.table();
@@ -2839 +2935 @@
     ArrayColumn<Float> specCol(tab, "SPECTRA");
     ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
+
     for (int i=0; i<int(tab.nrow()); ++i) {
       Vector<Float> spec = specCol(i);
       Vector<uChar> flag = flagCol(i);
-      int fstart = -1;
-      int fend = -1;
-      for (unsigned int k=0; k < flag.nelements(); ++k ) {
-        if (flag[k] > 0) {
-          fstart = k;
-          while (flag[k] > 0 && k < flag.nelements()) {
-            fend = k;
-            k++;
-          }
-        }
-        Float interp = 0.0;
-        if (fstart-1 > 0 ) {
-          interp = spec[fstart-1];
-          if (fend+1 < Int(spec.nelements())) {
-            interp = (interp+spec[fend+1])/2.0;
-          }
-        } else {
-          interp = spec[fend+1];
-        }
-        if (fstart > -1 && fend > -1) {
-          for (int j=fstart;j<=fend;++j) {
-            spec[j] = interp;
-          }
-        }
-        fstart =-1;
-        fend = -1;
-      }
+      doZeroOrderInterpolation(spec, flag);
+
       Vector<Complex> lags;
       ffts.fft0(lags, spec);
+
       Int lag0(start+0.5);
       Int lag1(end+0.5);
@@ -2878 +2951 @@
       }
       Int lstart =  max(0, lag0);
-      Int lend =  min(Int(lags.nelements()-1), lag1);
+      Int lend   =  min(Int(lags.nelements()-1), lag1);
       if (lstart == lend) {
         lags[lstart] = Complex(0.0);
@@ -2891 +2964 @@
         }
       }
+
       ffts.fft0(spec, lags);
+
       specCol.put(i, spec);
     }

trunk/src/STMath.h

@@ -13 +13 @@
 #define ASAPSTMATH_H
 
+#include <map>
 #include <string>
-#include <map>
 
 #include <casa/aips.h>
+#include <casa/Arrays/Vector.h>
+#include <casa/BasicSL/String.h>
 #include <casa/Utilities/CountedPtr.h>
-#include <casa/BasicSL/String.h>
-#include <casa/Arrays/Vector.h>
+
+#include <scimath/Mathematics/FFTServer.h>
 #include <scimath/Mathematics/InterpolateArray1D.h>
 
+#include "Logger.h"
 #include "Scantable.h"
 #include "STDefs.h"
 #include "STPol.h"
-#include "Logger.h"
 
 namespace asap {
@@ -345 +347 @@
     lagFlag( const casa::CountedPtr<Scantable>& in, double start,
              double end, const std::string& mode="frequency");
+
+  std::vector<float>
+    fft( const casa::CountedPtr<Scantable>& in,
+         const std::vector<int>& whichrow, 
+         bool getRealImag=false );
 
   // test for average spectra with different channel/resolution
@@ -426 +433 @@
                                         int index ) ;
   double getMJD( string strtime ) ;
+  void doFFT( std::vector<float>& out, 
+              casa::FFTServer< casa::Float, casa::Complex >& ffts, 
+              bool getRealImag, 
+              casa::Vector<casa::Float>& spec, 
+              casa::Vector<casa::uChar>& flag ) ;
+  void doZeroOrderInterpolation( casa::Vector<casa::Float>& spec, 
+                                 casa::Vector<casa::uChar>& flag) ;
 
   bool insitu_;

trunk/src/STMathWrapper.h

@@ -213 +213 @@
                                             mode)); }
 
+  std::vector<float> fft( const ScantableWrapper& in,
+                          const std::vector<int>& whichrow,
+                          const bool getRealImag=false )
+  { return STMath::fft(in.getCP(), whichrow, getRealImag); }
+
   // test for average spectra with different channel/resolution
   ScantableWrapper

trunk/src/python_STMath.cpp

@@ -77 +77 @@
         .def("_mx_extract", &STMathWrapper::mxExtract)
         .def("_lag_flag", &STMathWrapper::lagFlag)
+        .def("_fft", &STMathWrapper::fft)
         // testing average spectra with different channel/resolution
         .def("_new_average", &STMathWrapper::new_average)