X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=bfb55e8dc4695ee5fb451c2847c3185d5e63990a;hb=e4ffe82feebf1df2ac1dd14633818eb9a739863f;hp=c62eb24330d59e0c55645e6918977bc73f5d726a;hpb=2d39e823ba389fc68e5317c422b55be006094252;p=ctsim.git

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index c62eb24..bfb55e8 100644
--- a/libctsim/filter.cpp
+++ b/libctsim/filter.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: filter.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $
+**  $Id: filter.cpp,v 1.6 2000/07/02 18:21:39 kevin Exp $
 **
 **  This program is free software; you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License (version 2) as
@@ -37,42 +37,106 @@
  *   int filt_type	Type of filter wanted
  *   double bw		Bandwidth of filter
  *   double xmin, xmax	Filter limits
- *   int n		Number of points in filter
+ *   int n		Number of points in signal
  *   double param	General input parameter to filters
- *   int domain		FREQ or SPATIAL domain wanted
+ *   int domain		FREQUENCY or SPATIAL domain wanted
  *   int numint		Number if intervals for calculating discrete inverse fourier xform
  *			for spatial domain filters.  For ANALYTIC solutions, use numint = 0
  */
 
-SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numint)
+SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalLength, int n, double param, const char* domainName, int numIntegral = 0)
 {
+  m_vecFilter = NULL;
+  m_vecFourierCosTable = NULL;
+  m_vecFourierSinTable = NULL;
   m_idFilter = convertFilterNameToID (filterName);
+  if (m_idFilter == FILTER_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid Filter name ";
+    m_failMessage += filterName;
+    return;
+  }
+  m_idFilterMethod = convertFilterMethodNameToID (filterMethodName);
+  if (m_idFilterMethod == FILTER_METHOD_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid filter method name ";
+    m_failMessage += filterMethodName;
+    return;
+  }
   m_idDomain = convertDomainNameToID (domainName);
-  init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numint);
+  if (m_idDomain == DOMAIN_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid domain name ";
+    m_failMessage += domainName;
+    return;
+  }
+  init (m_idFilter, m_idFilterMethod, bw, signalLength, n, param, m_idDomain, numIntegral);
+}
+
+SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalLength, int n, double param, const DomainID domainID, int numIntegral = 0)
+{
+  init (filterID, filterMethodID, bw, signalLength, n, param, domainID, numIntegral);
 }
 
-SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
+SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0)
 {
-  init (filterID, bw, xmin, xmax, n, param, domainID, numint);
+  m_bw = bw;
+  m_nPoints = 0;
+  m_vecFilter = NULL;
+  m_vecFourierCosTable = NULL;
+  m_vecFourierSinTable = NULL;
+  m_filterParam = param;  
+  m_numIntegral = numIntegral;
+  m_idFilter = convertFilterNameToID (filterName);
+  if (m_idFilter == FILTER_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid Filter name ";
+    m_failMessage += filterName;
+    return;
+  }
+  m_idDomain = convertDomainNameToID (domainName);
+  if (m_idDomain == DOMAIN_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid domain name ";
+    m_failMessage += domainName;
+    return;
+  }
 }
 
 void
-SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
+SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalLength, int n, double param, const DomainID domainID, int numint)
 {
   m_bw = bw;
   m_idFilter = filterID;
   m_idDomain = domainID;
-  if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) {
+  m_idFilterMethod = filterMethodID;
+  if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID) {
     m_fail = true;
     return;
   }
   m_nameFilter = convertFilterIDToName (m_idFilter);
   m_nameDomain = convertDomainIDToName (m_idDomain);
+  m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
   m_fail = false;
-  m_nPoints = n;
-  m_xmin = xmin;
-  m_xmax = xmax;
+  m_nSignalPoints = n;
+  m_nFilterPoints = 2 * m_nSignalPoints - 1;
+
+  m_signalLength = signalLength;
+  m_xmin = -signalLength;
+  m_xmax = signalLength;
+  m_numIntegral = numint;
+  m_filterParam = param;  
   m_vecFilter = new double[n];
+  if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
+    int nFourier = n * n + 1;
+    double angleIncrement = (2. * PI) / n;
+    m_vecFourierCosTable = new double[ nFourier ];
+    m_vecFourierSinTable = new double[ nFourier ];
+    for (int i = 0; i < nFourier; i++) {
+      m_vecFourierCosTable[i] = cos (angleIncrement * i);
+      m_vecFourierSinTable[i] = sin (angleIncrement * i);
+    }
+  }
 
   double xinc = (m_xmax - m_xmin) / (m_nPoints - 1);
 
@@ -85,7 +149,7 @@ SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax
 
     for (int i = 1; i <= sidelen; i++ )
       m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
-  } else if (m_idDomain == DOMAIN_FREQ) {
+  } else if (m_idDomain == DOMAIN_FREQUENCY) {
     double x;
     int i;
     for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
@@ -99,7 +163,8 @@ SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax
       else
 	m_vecFilter[i] = spatialResponseCalc (x, param, numint);
   } else {
-      sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", m_idDomain);
+      m_failMessage = "Illegal domain name ";
+      m_failMessage += m_idDomain;
       m_fail = true;
   }
 }
@@ -107,13 +172,15 @@ SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax
 SignalFilter::~SignalFilter (void)
 {
     delete m_vecFilter;
+    delete m_vecFourierSinTable;
+    delete m_vecFourierCosTable;
 }
 
 
-SignalFilter::FilterID
+const SignalFilter::FilterID
 SignalFilter::convertFilterNameToID (const char *filterName)
 {
-  FilterID filterID;
+  FilterID filterID = FILTER_INVALID;
 
   if (strcasecmp (filterName, FILTER_BANDLIMIT_STR) == 0)
     filterID = FILTER_BANDLIMIT;
@@ -135,10 +202,6 @@ SignalFilter::convertFilterNameToID (const char *filterName)
     filterID = FILTER_ABS_COSINE;
   else if (strcasecmp (filterName, FILTER_SHEPP_STR) == 0)
     filterID = FILTER_SHEPP;
-  else {
-    sys_error(ERR_WARNING, "Invalid filter type %s\n", filterName);
-    filterID = FILTER_INVALID;
-  }
 
   return (filterID);
 }
@@ -172,17 +235,57 @@ SignalFilter::convertFilterIDToName (const FilterID filterID)
   return (name);
 }
       
+const SignalFilter::FilterMethodID
+SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName)
+{
+  FilterMethodID fmID = FILTER_METHOD_INVALID;
+
+  if (strcasecmp (filterMethodName, FILTER_METHOD_CONVOLUTION_STR) == 0)
+    fmID = FILTER_METHOD_CONVOLUTION;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FOURIER_STR) == 0)
+    fmID = FILTER_METHOD_FOURIER;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_STR) == 0)
+    fmID = FILTER_METHOD_FFT;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_ZEROPAD_2_STR) == 0)
+    fmID = FILTER_METHOD_FFT_ZEROPAD_2;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_ZEROPAD_4_STR) == 0)
+    fmID = FILTER_METHOD_FFT_ZEROPAD_4;
+  else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_ZEROPAD_6_STR) == 0)
+    fmID = FILTER_METHOD_FFT_ZEROPAD_6;
+
+  return (fmID);
+}
+
+const char *
+SignalFilter::convertFilterMethodIDToName (const FilterMethodID fmID)
+{
+  const char *name = "";
+
+  if (fmID == FILTER_METHOD_CONVOLUTION)
+    return (FILTER_METHOD_CONVOLUTION_STR);
+  else if (fmID == FILTER_METHOD_FOURIER)
+    return (FILTER_METHOD_FOURIER_STR);
+  else if (fmID == FILTER_METHOD_FFT)
+    return (FILTER_METHOD_FFT_STR);
+  else if (fmID == FILTER_METHOD_FFT_ZEROPAD_2)
+    return (FILTER_METHOD_FFT_ZEROPAD_2_STR);
+  else if (fmID == FILTER_METHOD_FFT_ZEROPAD_4)
+    return (FILTER_METHOD_FFT_ZEROPAD_4_STR);
+  else if (fmID == FILTER_METHOD_FFT_ZEROPAD_6)
+    return (FILTER_METHOD_FFT_ZEROPAD_6_STR);
+
+  return (name);
+}
+
 const SignalFilter::DomainID
 SignalFilter::convertDomainNameToID (const char* const domainName)
 {
-  DomainID dID;
+  DomainID dID = DOMAIN_INVALID;
 
   if (strcasecmp (domainName, DOMAIN_SPATIAL_STR) == 0)
     dID = DOMAIN_SPATIAL;
-  else if (strcasecmp (domainName, DOMAIN_FREQ_STR) == 0)
-    dID = DOMAIN_FREQ;
-  else
-    dID = DOMAIN_INVALID;
+  else if (strcasecmp (domainName, DOMAIN_FREQUENCY_STR) == 0)
+    dID = DOMAIN_FREQUENCY;
 
   return (dID);
 }
@@ -194,28 +297,62 @@ SignalFilter::convertDomainIDToName (const DomainID domain)
 
   if (domain == DOMAIN_SPATIAL)
     return (DOMAIN_SPATIAL_STR);
-  else if (domain == DOMAIN_FREQ)
-    return (DOMAIN_FREQ_STR);
+  else if (domain == DOMAIN_FREQUENCY)
+    return (DOMAIN_FREQUENCY_STR);
 
   return (name);
 }
 
 
+void
+SignalFilter::filterSignal (const double input[], double output[], double dx, const int n) const
+{
+  if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
+    for (int i = 0; i < n; i++)
+      output[i] = convolve (input, dx, i, n);
+  } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
+    complex<double> fftSignal[n];
+    complex<double> complexOutput;
+    finiteFourierTransform (input, fftSignal, 1);
+    finiteFourierTransform (fftSignal, complexOutput, -1);
+    for (int i = 0; i < n; i++)
+      output[i] = complexOutput[i].mag();
+  }
+}
+
+void
+SignalFilter::filterSignal (const float input[], double output[], double dx, const int n) const
+{
+  if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
+    for (int i = 0; i < n; i++)
+      output[i] = convolve (input, dx, i, n);
+  }
+}
+
+
 double
-SignalFilter::response (const char* filterName, const char* domainName, double bw, double x, double filt_param)
+SignalFilter::response (double x)
 {
   double response = 0;
-  FilterID filterID = convertFilterNameToID (filterName);
-  DomainID domainID = convertDomainNameToID (domainName);
 
-  if (domainID == DOMAIN_SPATIAL)
-    response = spatialResponseAnalytic (filterID, bw, x, filt_param);
-  else if (domainID == DOMAIN_FREQ)
-    response = frequencyResponse (filterID, bw, x, filt_param);
+  if (m_idDomain == DOMAIN_SPATIAL)
+    response = spatialResponse (m_idFilter, m_bw, x, m_filterParam, m_numIntegral);
+  else if (m_idDomain == DOMAIN_FREQUENCY)
+    response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam);
 
   return (response);
 }
 
+
+double 
+SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double param, int nIntegral = 0)
+{
+  if (nIntegral == 0)
+    return spatialResponseAnalytic (filterID, bw, x, param);
+  else
+    return spatialResponseCalc (filterID, bw, x, param, nIntegral);
+}
+
 /* NAME
  *   filter_spatial_response_calc	Calculate filter by discrete inverse fourier
  *				        transform of filters's frequency
@@ -232,9 +369,9 @@ SignalFilter::response (const char* filterName, const char* domainName, double b
  */
 
 double 
-SignalFilter::spatialResponseCalc (double x, double param, int n) const
+SignalFilter::spatialResponseCalc (double x, double param, int nIntegral) const
 {
-  return (spatialResponseCalc (m_idFilter, m_bw, x, param, n));
+  return (spatialResponseCalc (m_idFilter, m_bw, x, param, nIntegral));
 }
 
 double 
@@ -451,10 +588,9 @@ SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, d
 double 
 SignalFilter::integral_abscos (double u, double w)
 {
-  if (fabs (u) > F_EPSILON)
-    return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
-  else
-    return (w * w / 2);
+  return (fabs (u) > F_EPSILON 
+     ? (cos(u * w) - 1) / (u * u) + w / u * sin (u * w) 
+     : (w * w / 2));
 }
 
 
@@ -519,3 +655,57 @@ for (int i = 0; i < np; i++)
   return (sum * dx);
 }
 
+
+void
+SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
+{
+  if (direction < 0)
+    direction = -1;
+  else 
+    direction = 1;
+    
+  double angleIncrement = 2 * PI / n;
+  for (int i = 0; i < n; i++) {
+    double sumReal = 0;
+    double sumImag = 0;
+    for (int j = 0; j < n; j++) {
+      double angle = i * j * angleIncrement * direction;
+      sumReal += input[i] * cos(angle);
+      sumImag += input[i] * sin(angle);
+    }
+    if (direction > 0) {
+      sumReal /= n;
+      sumImag /= n;
+    }
+    output[i] = complex<double> (sumReal, sumImag);
+  }
+}
+
+void
+SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
+{
+  if (direction < 0)
+    direction = -1;
+  else 
+    direction = 1;
+    
+  double angleIncrement = 2 * PI / m_nPoints;
+  for (int i = 0; i < m_nPoints; i++) {
+    double sumReal = 0, sumImag = 0;
+    for (int j = 0; j < m_nPoints; j++) {
+      int tableIndex = i * j;
+      if (direction > 0) {
+	sumReal += input[i] * m_vecFourierCosTable[tableIndex];
+	sumImag += input[i] * m_vecFourierSinTable[tableIndex];
+      } else {
+	sumReal += input[i] * m_vecFourierCosTable[tableIndex];
+	sumImag -= input[i] * m_vecFourierSinTable[tableIndex];
+      }
+    }
+    if (direction > 0) {
+      sumReal /= m_nPoints;
+      sumImag /= m_nPoints;
+    }
+    output[i] = complex<double> (sumReal, sumImag);
+  }
+}