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

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index b2d3062..c62eb24 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.2 2000/06/20 17:54:51 kevin Exp $
+**  $Id: filter.cpp,v 1.3 2000/06/22 10:17:28 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
@@ -29,63 +29,192 @@
 
 
 /* NAME
- *   filter_generate				Generate a filter
+ *   SignalFilter::SignalFilter     Construct a signal
  *
  * SYNOPSIS
- *   f = filter_generate (filt_type, bw, xmin, xmax, n, param, domain, analytic)
- *   double f				Generated filter vector
- *   int filt_type			Type of filter wanted
- *   double bw				Bandwidth of filter
- *   double xmin, xmax			Filter limits
- *   int n				Number of points in filter
- *   double param			General input parameter to filters
- *   int domain				FREQ 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
+ *   f = SignalFilter (filt_type, bw, xmin, xmax, n, param, domain, analytic)
+ *   double f		Generated filter vector
+ *   int filt_type	Type of filter wanted
+ *   double bw		Bandwidth of filter
+ *   double xmin, xmax	Filter limits
+ *   int n		Number of points in filter
+ *   double param	General input parameter to filters
+ *   int domain		FREQ 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 FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint)
+SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numint)
+{
+  m_idFilter = convertFilterNameToID (filterName);
+  m_idDomain = convertDomainNameToID (domainName);
+  init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numint);
+}
+
+SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
+{
+  init (filterID, bw, xmin, xmax, n, param, domainID, numint);
+}
+
+void
+SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
 {
-  m_vecFilter = new double [n];
-  m_filterType = filt_type;
   m_bw = bw;
+  m_idFilter = filterID;
+  m_idDomain = domainID;
+  if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) {
+    m_fail = true;
+    return;
+  }
+  m_nameFilter = convertFilterIDToName (m_idFilter);
+  m_nameDomain = convertDomainIDToName (m_idDomain);
+  m_fail = false;
+  m_nPoints = n;
+  m_xmin = xmin;
+  m_xmax = xmax;
+  m_vecFilter = new double[n];
 
-  double xinc = (xmax - xmin) / (n - 1);
+  double xinc = (m_xmax - m_xmin) / (m_nPoints - 1);
 
-  if (m_filterType == FILTER_SHEPP) {
+  if (m_idFilter == FILTER_SHEPP) {
     double a = 2 * m_bw;
     double c = - 4. / (a * a);
-    int center = (n - 1) / 2;
+    int center = (m_nPoints - 1) / 2;
     int sidelen = center;
     m_vecFilter[center] = 4. / (a * a);
-    
+
     for (int i = 1; i <= sidelen; i++ )
       m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
-  } else if (domain == D_FREQ) {
+  } else if (m_idDomain == DOMAIN_FREQ) {
     double x;
     int i;
-    for (x = xmin, i = 0; i < n; x += xinc, i++)
+    for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
       m_vecFilter[i] = frequencyResponse (x, param);
-  } else if (domain == D_SPATIAL) {
+  } else if (m_idDomain == DOMAIN_SPATIAL) {
     double x;
     int i;
-    for (x = xmin, i = 0; i < n; x += xinc, i++)
+    for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
       if (numint == 0)
 	m_vecFilter[i] = spatialResponseAnalytic (x, param);
       else
 	m_vecFilter[i] = spatialResponseCalc (x, param, numint);
   } else {
-    sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain);
+      sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", m_idDomain);
+      m_fail = true;
   }
 }
 
 SignalFilter::~SignalFilter (void)
 {
-  delete m_vecFilter;
+    delete m_vecFilter;
+}
+
+
+SignalFilter::FilterID
+SignalFilter::convertFilterNameToID (const char *filterName)
+{
+  FilterID filterID;
+
+  if (strcasecmp (filterName, FILTER_BANDLIMIT_STR) == 0)
+    filterID = FILTER_BANDLIMIT;
+  else if (strcasecmp (filterName, FILTER_HAMMING_STR) == 0)
+    filterID = FILTER_G_HAMMING;
+  else if (strcasecmp (filterName, FILTER_SINC_STR) == 0)
+    filterID = FILTER_SINC;
+  else if (strcasecmp (filterName, FILTER_COS_STR) == 0)
+    filterID = FILTER_COSINE;
+  else if (strcasecmp (filterName, FILTER_TRIANGLE_STR) == 0)
+    filterID = FILTER_TRIANGLE;
+  else if (strcasecmp (filterName, FILTER_ABS_BANDLIMIT_STR) == 0)
+    filterID = FILTER_ABS_BANDLIMIT;
+  else if (strcasecmp (filterName, FILTER_ABS_HAMMING_STR) == 0)
+    filterID = FILTER_ABS_G_HAMMING;
+  else if (strcasecmp (filterName, FILTER_ABS_SINC_STR) == 0)
+    filterID = FILTER_ABS_SINC;
+  else if (strcasecmp (filterName, FILTER_ABS_COS_STR) == 0)
+    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);
+}
+
+const char *
+SignalFilter::convertFilterIDToName (const FilterID filterID)
+{
+  const char *name = "";
+
+  if (filterID == FILTER_SHEPP)
+    name = FILTER_SHEPP_STR;
+  else if (filterID == FILTER_ABS_COSINE)
+    name = FILTER_ABS_COS_STR;
+  else if (filterID == FILTER_ABS_SINC)
+    name = FILTER_ABS_SINC_STR;
+  else if (filterID == FILTER_ABS_G_HAMMING)
+    name = FILTER_ABS_HAMMING_STR;
+  else if (filterID == FILTER_ABS_BANDLIMIT)
+    name = FILTER_ABS_BANDLIMIT_STR;
+  else if (filterID == FILTER_COSINE)
+    name = FILTER_COS_STR;
+  else if (filterID == FILTER_SINC)
+    name = FILTER_SINC_STR;
+  else if (filterID == FILTER_G_HAMMING)
+    name = FILTER_HAMMING_STR;
+  else if (filterID == FILTER_BANDLIMIT)
+    name = FILTER_BANDLIMIT_STR;
+  else if (filterID == FILTER_TRIANGLE)
+    name = FILTER_TRIANGLE_STR;
+	    
+  return (name);
+}
+      
+const SignalFilter::DomainID
+SignalFilter::convertDomainNameToID (const char* const domainName)
+{
+  DomainID dID;
+
+  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;
+
+  return (dID);
 }
 
+const char *
+SignalFilter::convertDomainIDToName (const DomainID domain)
+{
+  const char *name = "";
+
+  if (domain == DOMAIN_SPATIAL)
+    return (DOMAIN_SPATIAL_STR);
+  else if (domain == DOMAIN_FREQ)
+    return (DOMAIN_FREQ_STR);
+
+  return (name);
+}
+
+
+double
+SignalFilter::response (const char* filterName, const char* domainName, double bw, double x, double filt_param)
+{
+  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);
+
+  return (response);
+}
 
 /* NAME
  *   filter_spatial_response_calc	Calculate filter by discrete inverse fourier
@@ -105,15 +234,15 @@ SignalFilter::~SignalFilter (void)
 double 
 SignalFilter::spatialResponseCalc (double x, double param, int n) const
 {
-  return (spatialResponseCalc (m_filterType, m_bw, x, param, n));
+  return (spatialResponseCalc (m_idFilter, m_bw, x, param, n));
 }
 
 double 
-SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double param, int n)
+SignalFilter::spatialResponseCalc (FilterID filterID, double bw, double x, double param, int n)
 {
   double zmin, zmax;
 
-  if (fType == FILTER_TRIANGLE) {
+  if (filterID == FILTER_TRIANGLE) {
     zmin = 0;
     zmax = bw;
   } else {
@@ -125,7 +254,7 @@ SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double
   double z = zmin;
   double q [n];
   for (int i = 0; i < n; i++, z += zinc)
-    q[i] = frequencyResponse (fType, bw, z, param) * cos (TWOPI * z * x);
+    q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
   
   double y = 2 * integrateSimpson (zmin, zmax, q, n);
   
@@ -148,17 +277,17 @@ SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double
 double 
 SignalFilter::frequencyResponse (double u, double param) const
 {
-  return frequencyResponse (m_filterType, m_bw, u, param);
+  return frequencyResponse (m_idFilter, m_bw, u, param);
 }
 
 
 double 
-SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double param)
+SignalFilter::frequencyResponse (FilterID filterID, double bw, double u, double param)
 {
   double q;
   double au = fabs (u);
 
-  switch (fType) {
+  switch (filterID) {
   case FILTER_BANDLIMIT:
     if (au >= bw / 2)
       q = 0.;
@@ -209,7 +338,7 @@ SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double p
     break;
   default:
     q = 0;
-    sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", fType);
+    sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
     break;
   }
   return (q);
@@ -234,11 +363,11 @@ SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double p
 double 
 SignalFilter::spatialResponseAnalytic (double x, double param) const
 {
-  return spatialResponseAnalytic (m_filterType, m_bw, x, param);
+  return spatialResponseAnalytic (m_idFilter, m_bw, x, param);
 }
 
 double 
-SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, double param)
+SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, double param)
 {
   double q, temp;
   double u = TWOPI * x;
@@ -246,7 +375,7 @@ SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, do
   double b = PI / bw;
   double b2 = TWOPI / bw;
 
-  switch (fType) {
+  switch (filterID) {
   case FILTER_BANDLIMIT:
     q = bw * sinc(u * w, 1.0);
     break;
@@ -284,7 +413,7 @@ SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, do
     break;
   case FILTER_ABS_SINC:
   default:
-    sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", fType);
+    sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
     q = 0;
     break;
   }
@@ -356,68 +485,36 @@ SignalFilter::integral_abscos (double u, double w)
  */
 
 double 
-SignalFilter::convolve (const double func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
+SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const
 {
   double sum = 0.0;
 
-  if (func_type == FUNC_BOTH) {
 #if UNOPTIMIZED_CONVOLUTION
-    for (int i = 0; i < np; i++)
-      sum += func[i] * m_vecFilter[n - i + (np - 1)];
+  for (int i = 0; i < np; i++)
+    sum += func[i] * m_vecFilter[n - i + (np - 1)];
 #else
-    double* f2 = m_vecFilter + n + (np - 1);
-    for (int i = 0; i < np; i++)
-      sum += *func++ * *f2--;
+  double* f2 = m_vecFilter + n + (np - 1);
+  for (int i = 0; i < np; i++)
+    sum += *func++ * *f2--;
 #endif
-  }  else if (func_type == FUNC_EVEN) {
-    for (int i = 0; i < np; i++) {
-      int k = abs (n - i);
-      sum += func[i] * m_vecFilter[k];
-    }
-  } else if (func_type == FUNC_ODD) {
-    for (int i = 0; i < np; i++) {
-      int k = n - i;
-      if (k < 0)
-	sum -= func[i] * m_vecFilter[k];
-      else
-	sum += func[i] * m_vecFilter[k];
-    }
-  } else
-    sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
 
   return (sum * dx);
 }
 
 
 double 
-SignalFilter::convolve (const float func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
+SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const
 {
   double sum = 0.0;
 
-  if (func_type == FUNC_BOTH) {
 #if UNOPTIMIZED_CONVOLUTION
-    for (int i = 0; i < np; i++)
-      sum += func[i] * m_vecFilter[n - i + (np - 1)];
+for (int i = 0; i < np; i++)
+  sum += func[i] * m_vecFilter[n - i + (np - 1)];
 #else
-    double* f2 = m_vecFilter + n + (np - 1);
-    for (int i = 0; i < np; i++)
-      sum += *func++ * *f2--;
+double* f2 = m_vecFilter + n + (np - 1);
+for (int i = 0; i < np; i++)
+  sum += *func++ * *f2--;
 #endif
-  }  else if (func_type == FUNC_EVEN) {
-    for (int i = 0; i < np; i++) {
-      int k = abs (n - i);
-      sum += func[i] * m_vecFilter[k];
-    }
-  } else if (func_type == FUNC_ODD) {
-    for (int i = 0; i < np; i++) {
-      int k = n - i;
-      if (k < 0)
-	sum -= func[i] * m_vecFilter[k];
-      else
-	sum += func[i] * m_vecFilter[k];
-    }
-  } else
-    sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
 
   return (sum * dx);
 }