include/filter.h

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **      Name:         filter.h
   5 **      Purpose:      Signal filter header file
   6 **      Programmer:   Kevin Rosenberg
   7 **      Date Started: June 2000
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: filter.h,v 1.18 2000/08/03 09:57:33 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #ifndef FILTER_H
  29 #define FILTER_H
  30
  31
  32 #ifdef HAVE_CONFIG_H
  33 #include "config.h"
  34 #endif
  35 #ifdef HAVE_FFTW
  36 #include <fftw.h>
  37 #include <rfftw.h>
  38 #endif
  39
  40 #include <complex>
  41
  42
  43 class SignalFilter {
  44  public:
  45
  46     static const int FILTER_INVALID;
  47     static const int FILTER_ABS_BANDLIMIT;      // filter times |x|
  48     static const int FILTER_ABS_SINC;
  49     static const int FILTER_ABS_G_HAMMING;
  50     static const int FILTER_ABS_COSINE;
  51     static const int FILTER_SHEPP;
  52     static const int FILTER_BANDLIMIT;
  53     static const int FILTER_SINC;
  54     static const int FILTER_G_HAMMING;
  55     static const int FILTER_COSINE;
  56     static const int FILTER_TRIANGLE;
  57
  58     static const int FILTER_METHOD_INVALID;
  59     static const int FILTER_METHOD_CONVOLUTION;
  60     static const int FILTER_METHOD_FOURIER;
  61     static const int FILTER_METHOD_FOURIER_TABLE;
  62     static const int FILTER_METHOD_FFT;
  63 #if HAVE_FFTW
  64     static const int FILTER_METHOD_FFTW;
  65     static const int FILTER_METHOD_RFFTW;
  66 #endif
  67
  68     static const int DOMAIN_INVALID;
  69     static const int DOMAIN_FREQUENCY;
  70     static const int DOMAIN_SPATIAL;
  71
  72     static const int FREQUENCY_FILTER_INVALID;
  73     static const int FREQUENCY_FILTER_DIRECT_FREQUENCY;
  74     static const int FREQUENCY_FILTER_INVERSE_SPATIAL;
  75
  76     SignalFilter (const char* filterName, const char* filterMethodName,double bw, double signalIncrement, int n, double param, const char* domainName, const char* frequencyFilterName, const int zeropad = 0, const int preinterpolationFactor = 1);
  77
  78     SignalFilter (const int filt_type, int filterMethodID, double bw, double signalIncrement, int n, double param, const int domain, int filterFilterID, const int zeropad = 0, const int preinterpolationFactor = 1);
  79
  80     SignalFilter (const char* filterName, const char* domainName, double bw, double param);
  81
  82     ~SignalFilter (void);
  83
  84     double* getFilter (void) const
  85       { return m_vecFilter; }
  86
  87     int getNFilterPoints (void) const
  88         { return m_nFilterPoints; }
  89
  90     double convolve (const double f[], const double dx, const int n, const int np) const;
  91
  92     double convolve (const float f[], const double dx, const int n, const int np) const;
  93
  94     void filterSignal (const double input[], double output[]) const;
  95     void filterSignal (const float input[], double output[]) const;
  96
  97     static void finiteFourierTransform (const double input[], complex<double> output[], const int n, const int direction);
  98     static void finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, const int direction);
  99     static void finiteFourierTransform (const complex<double> input[], double output[], const int n, const int direction);
 100
 101     void finiteFourierTransform (const double input[], complex<double> output[], const int direction) const;
 102     void finiteFourierTransform (const complex<double> input[], complex<double> output[], const int direction) const;
 103     void finiteFourierTransform (const complex<double> input[], double output[], const int direction) const;
 104
 105     void setTraceLevel (int traceLevel) {m_traceLevel = traceLevel; }
 106
 107     bool fail(void) const       {return m_fail;}
 108     const string& failMessage(void) const {return m_failMessage;}
 109
 110     const string& nameFilter(void) const        { return m_nameFilter;}
 111     const string& nameDomain(void) const        { return m_nameDomain;}
 112     const int idFilter(void) const      { return m_idFilter;}
 113     const int idDomain(void) const      { return m_idDomain;}
 114     const int idFrequencyFilter() const { return m_idFrequencyFilter;}
 115     const double getFilterMin(void) const {return m_filterMin;}
 116     const double getFilterMax(void) const {return m_filterMax;}
 117     const double getFilterIncrement(void) const {return m_filterInc;}
 118
 119     double response (double x);
 120
 121     static double spatialResponse (int fType, double bw, double x, double param);
 122
 123     static double frequencyResponse (int fType, double bw, double u, double param);
 124
 125     static double spatialResponseAnalytic (int fType, double bw, double x, double param);
 126
 127     static double spatialResponseCalc (int fType, double bw, double x, double param, int nIntegral);
 128
 129     static void setNumIntegral(int nIntegral) {N_INTEGRAL = nIntegral;}
 130
 131   static const int getFilterCount() {return s_iFilterCount;}
 132   static const char** getFilterNameArray() {return s_aszFilterName;}
 133   static const char** getFilterTitleArray() {return s_aszFilterTitle;}
 134   static int convertFilterNameToID (const char* const filterName);
 135   static const char* convertFilterIDToName (const int idFilter);
 136   static const char* convertFilterIDToTitle (const int idFilter);
 137
 138   static const int getFilterMethodCount() {return s_iFilterMethodCount;}
 139   static const char** getFilterMethodNameArray() {return s_aszFilterMethodName;}
 140   static const char** getFilterMethodTitleArray() {return s_aszFilterMethodTitle;}
 141   static int convertFilterMethodNameToID (const char* const filterMethodName);
 142   static const char* convertFilterMethodIDToName (const int idFilterMethod);
 143   static const char* convertFilterMethodIDToTitle (const int idFilterMethod);
 144
 145   static const int getDomainCount() {return s_iDomainCount;}
 146   static const char** getDomainNameArray() {return s_aszDomainName;}
 147   static const char** getDomainTitleArray() {return s_aszDomainTitle;}
 148   static int convertDomainNameToID (const char* const domainName);
 149   static const char* convertDomainIDToName (const int idDomain);
 150   static const char* convertDomainIDToTitle (const int idDomain);
 151   static int convertFrequencyFilterNameToID (const char* const ffName);
 152   static const char* convertFrequencyFilterIDToName (const int idFF);
 153   static const char* convertFrequencyFilterIDToTitle (const int idFF);
 154
 155  private:
 156     double m_bw;
 157     int m_nFilterPoints;
 158     int m_nSignalPoints;
 159     double m_signalInc;
 160     double m_filterMin;
 161     double m_filterMax;
 162     double m_filterInc;
 163     double* m_vecFilter;
 164     double* m_vecFourierCosTable;
 165     double* m_vecFourierSinTable;
 166     complex<double>* m_complexVecFilter;
 167 #ifdef HAVE_FFTW
 168     fftw_real* m_vecRealFftInput, *m_vecRealFftSignal;
 169     rfftw_plan m_realPlanForward, m_realPlanBackward;
 170     fftw_complex* m_vecComplexFftInput, *m_vecComplexFftSignal;
 171     fftw_plan m_complexPlanForward, m_complexPlanBackward;
 172 #endif
 173
 174     bool m_fail;
 175     string m_failMessage;
 176     string m_nameFilter;
 177     string m_nameFilterMethod;
 178     string m_nameDomain;
 179     string m_nameFrequencyFilter;
 180     int m_idFilter;
 181     int m_idFilterMethod;
 182     int m_idFrequencyFilter;
 183     int m_idDomain;
 184     double m_filterParam;
 185     int m_traceLevel;
 186     int m_zeropad;
 187     int m_nOutputPoints;
 188     int m_preinterpolationFactor;
 189
 190     static const char* s_aszFilterName[];
 191     static const char* s_aszFilterTitle[];
 192     static const int s_iFilterCount;
 193     static const char* s_aszFilterMethodName[];
 194     static const char* s_aszFilterMethodTitle[];
 195     static const int s_iFilterMethodCount;
 196     static const char* s_aszDomainName[];
 197     static const char* s_aszDomainTitle[];
 198     static const int s_iDomainCount;
 199     static const char* s_aszFrequencyFilterName[];
 200     static const char* s_aszFrequencyFilterTitle[];
 201     static const int s_iFrequencyFilterCount;
 202
 203     static int N_INTEGRAL;
 204
 205     static const bool haveAnalyticSpatial (const int filterID);
 206
 207     void init (const int filt_type, const int filterMethod, double bw, double signalIncrement, int n, double param, const int domain, const int frequencyFilter, const int zeropad, const int preInterpScale);
 208
 209     double spatialResponseCalc (double x, double param) const;
 210
 211     double spatialResponseAnalytic (double x, double param) const;
 212
 213     double frequencyResponse (double u, double param) const;
 214
 215     static double sinc (double x, double mult)
 216       { return (fabs(x) > F_EPSILON ? (sin (x * mult) / x) : 1.0); }
 217
 218     static double integral_abscos (double u, double w);
 219
 220 };
 221
 222 #endif