X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=d7d8b973c61a41f233d5a03a7c1dd3b894a5f8fc;hb=d158a6da6ec776fe98961f73f2ea74cf55de4700;hp=1eda2366ef5da6f70616562edb32a15a2df0c1f6;hpb=286d655a25df2668bd65ad365676c6ecc94415a1;p=ctsim.git

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 1eda236..d7d8b97 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.20 2000/07/22 15:45:33 kevin Exp $
+**  $Id: filter.cpp,v 1.24 2000/08/03 09:57:33 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
@@ -30,41 +30,41 @@
 int SignalFilter::N_INTEGRAL=500;  //static member
 
 const int SignalFilter::FILTER_INVALID = -1 ;
-const int SignalFilter::FILTER_BANDLIMIT = 0;
-const int SignalFilter::FILTER_SINC = 1;
-const int SignalFilter::FILTER_G_HAMMING = 2;
-const int SignalFilter::FILTER_COSINE = 3;
-const int SignalFilter::FILTER_TRIANGLE = 4;
-const int SignalFilter::FILTER_ABS_BANDLIMIT = 5; 	// filter times |x = |
-const int SignalFilter::FILTER_ABS_SINC = 6;
-const int SignalFilter::FILTER_ABS_G_HAMMING = 7;
-const int SignalFilter::FILTER_ABS_COSINE = 8;
-const int SignalFilter::FILTER_SHEPP = 9;
+const int SignalFilter::FILTER_ABS_BANDLIMIT = 0; 	// filter times |x = |
+const int SignalFilter::FILTER_ABS_SINC = 1;
+const int SignalFilter::FILTER_ABS_G_HAMMING = 2;
+const int SignalFilter::FILTER_ABS_COSINE = 3;
+const int SignalFilter::FILTER_SHEPP = 4;
+const int SignalFilter::FILTER_BANDLIMIT = 5;
+const int SignalFilter::FILTER_SINC = 6;
+const int SignalFilter::FILTER_G_HAMMING = 7;
+const int SignalFilter::FILTER_COSINE = 8;
+const int SignalFilter::FILTER_TRIANGLE = 9;
 
 const char* SignalFilter::s_aszFilterName[] = {
-  {"bandlimit"},
-  {"sinc"},
-  {"hamming"},
-  {"cosine"},
-  {"triangle"},
   {"abs_bandlimit"},
   {"abs_sinc"},
   {"abs_hamming"},
   {"abs_cosine"},
   {"shepp"},
+  {"bandlimit"},
+  {"sinc"},
+  {"hamming"},
+  {"cosine"},
+  {"triangle"},
 };
 
 const char* SignalFilter::s_aszFilterTitle[] = {
-  {"Bandlimit"},
-  {"Sinc"},
-  {"Hamming"},
-  {"Cosine"},
-  {"Triangle"},
   {"Abs(w) * Bandlimit"},
   {"Abs(w) * Sinc"},
   {"Abs(w) * Hamming"},
   {"Abs(w) * Cosine"},
   {"Shepp"},
+  {"Bandlimit"},
+  {"Sinc"},
+  {"Hamming"},
+  {"Cosine"},
+  {"Triangle"},
 };
 
 const int SignalFilter::s_iFilterCount = sizeof(s_aszFilterName) / sizeof(const char*);
@@ -121,6 +121,23 @@ const char* SignalFilter::s_aszDomainTitle[] = {
 const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const char*);
 
 
+const int SignalFilter::FREQUENCY_FILTER_INVALID = -1;
+const int SignalFilter::FREQUENCY_FILTER_DIRECT_FREQUENCY = 0;
+const int SignalFilter::FREQUENCY_FILTER_INVERSE_SPATIAL = 1;
+    
+const char* SignalFilter::s_aszFrequencyFilterName[] = {
+  {"direct_frequency"},
+  {"inverse_spatial"},
+};
+
+const char* SignalFilter::s_aszFrequencyFilterTitle[] = {
+  {"Direct Frequency"},
+  {"Inverse Spatial"},
+};
+
+const int SignalFilter::s_iFrequencyFilterCount = sizeof(s_aszFrequencyFilterName) / sizeof(const char*);
+
+
 /* NAME
  *   SignalFilter::SignalFilter     Construct a signal
  *
@@ -135,11 +152,9 @@ const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const
  *   int domain		FREQUENCY or SPATIAL domain wanted
  */
 
-SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, int zeropad = 0, int preinterpolationFactor = 1)
+SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, const char* frequencyFilterName, int zeropad = 0, int preinterpolationFactor = 1)
+  : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false)
 {
-  m_vecFilter = NULL;
-  m_vecFourierCosTable = NULL;
-  m_vecFourierSinTable = NULL;
   m_idFilter = convertFilterNameToID (filterName);
   if (m_idFilter == FILTER_INVALID) {
     m_fail = true;
@@ -161,22 +176,28 @@ SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName
     m_failMessage += domainName;
     return;
   }
-  init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, zeropad, preinterpolationFactor);
+  m_idFrequencyFilter = convertFrequencyFilterNameToID (frequencyFilterName);
+  if (m_idFrequencyFilter == FREQUENCY_FILTER_INVALID) {
+    m_fail = true;
+    m_failMessage = "Invalid frequency filter name ";
+    m_failMessage += frequencyFilterName;
+    return;
+  }
+  init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, m_idFrequencyFilter, zeropad, preinterpolationFactor);
 }
 
-SignalFilter::SignalFilter (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const int domainID, int zeropad = 0, int preinterpolationFactor = 1)
+SignalFilter::SignalFilter (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const int domainID, int frequencyFilterID, int zeropad = 0, int preinterpolationFactor = 1)
+  : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false)
 {
-  init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, zeropad, preinterpolationFactor);
+  init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, frequencyFilterID, zeropad, preinterpolationFactor);
 }
 
 SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param)
+  : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false)
 {
   m_bw = bw;
   m_nSignalPoints = 0;
   m_nFilterPoints = 0;
-  m_vecFilter = NULL;
-  m_vecFourierCosTable = NULL;
-  m_vecFourierSinTable = NULL;
   m_filterParam = param;  
   m_idFilter = convertFilterNameToID (filterName);
   if (m_idFilter == FILTER_INVALID) {
@@ -195,13 +216,14 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub
 }
 
 void
-SignalFilter::init (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double filterParam, const int domainID, int zeropad, int preinterpolationFactor)
+SignalFilter::init (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double filterParam, const int domainID, const int frequencyFilterID, int zeropad, int preinterpolationFactor)
 {
   m_bw = bw;
   m_idFilter = filterID;
   m_idDomain = domainID;
   m_idFilterMethod = filterMethodID;
-  if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID) {
+  m_idFrequencyFilter = frequencyFilterID;
+  if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFrequencyFilter == FREQUENCY_FILTER_INVALID) {
     m_fail = true;
     return;
   }
@@ -209,7 +231,7 @@ SignalFilter::init (const int filterID, const int filterMethodID, double bw, dou
   m_nameFilter = convertFilterIDToName (m_idFilter);
   m_nameDomain = convertDomainIDToName (m_idDomain);
   m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
-  m_fail = false;
+  m_nameFrequencyFilter = convertFrequencyFilterIDToName (m_idFrequencyFilter);
   m_nSignalPoints = nSignalPoints;
   m_signalInc = signalIncrement;
   m_filterParam = filterParam;  
@@ -321,11 +343,16 @@ SignalFilter::init (const int filterID, const int filterMethodID, double bw, dou
     } else if (m_idDomain == DOMAIN_SPATIAL) {
       double x;
       int i;
-      for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++)
+      for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) {
 	if (haveAnalyticSpatial(m_idFilter))
 	  m_vecFilter[i] = spatialResponseAnalytic (x, m_filterParam);
 	else
 	  m_vecFilter[i] = spatialResponseCalc (x, m_filterParam);
+#if LIMIT_BANDWIDTH_TRIAL
+	if (i < m_nFilterPoints / 4 || i > (m_nFilterPoints * 3) / 4)
+	  m_vecFilter[i] = 0;
+#endif
+      }
     } else {
       m_failMessage = "Illegal domain name ";
       m_failMessage += m_idDomain;
@@ -363,10 +390,10 @@ SignalFilter::convertFilterNameToID (const char *filterName)
   int filterID = FILTER_INVALID;
 
   for (int i = 0; i < s_iFilterCount; i++)
-      if (strcasecmp (filterName, s_aszFilterName[i]) == 0) {
-	filterID = i;
-	break;
-      }
+    if (strcasecmp (filterName, s_aszFilterName[i]) == 0) {
+      filterID = i;
+      break;
+    }
 
   return (filterID);
 }
@@ -413,7 +440,7 @@ SignalFilter::convertFilterMethodIDToName (const int fmID)
   static const char *name = "";
 
   if (fmID >= 0 && fmID < s_iFilterMethodCount)
-      return (s_aszFilterName [fmID]);
+      return (s_aszFilterMethodName [fmID]);
 
   return (name);
 }
@@ -465,6 +492,42 @@ SignalFilter::convertDomainIDToTitle (const int domainID)
   return (title);
 }
 
+int
+SignalFilter::convertFrequencyFilterNameToID (const char* const ffName)
+{
+  int ffID = FREQUENCY_FILTER_INVALID;
+
+  for (int i = 0; i < s_iFrequencyFilterCount; i++)
+   if (strcasecmp (ffName, s_aszFrequencyFilterName[i]) == 0) {
+      ffID = i;
+      break;
+    }
+
+  return (ffID);
+}
+
+const char *
+SignalFilter::convertFrequencyFilterIDToName (const int ffID)
+{
+  static const char *name = "";
+
+  if (ffID >= 0 && ffID < s_iFrequencyFilterCount)
+      return (s_aszFrequencyFilterName [ffID]);
+
+  return (name);
+}
+
+const char *
+SignalFilter::convertFrequencyFilterIDToTitle (const int ffID)
+{
+  static const char *name = "";
+
+  if (ffID >= 0 && ffID < s_iFrequencyFilterCount)
+      return (s_aszFrequencyFilterTitle [ffID]);
+
+  return (name);
+}
+
 void
 SignalFilter::filterSignal (const float input[], double output[]) const
 {