X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=221f617f869a716f759ef49349db6253442e5041;hp=6e5910a990030451664ba29a9680069488458628;hb=1a050c98763fbbc0662731b0b76953acede6f5d7;hpb=c8b19dfaffba9f06d8b6c40cb1bb83a8964867f7 diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index 6e5910a..221f617 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -1,9 +1,9 @@ /***************************************************************************** ** File IDENTIFICATION -** +** ** Name: filter.cpp ** Purpose: Routines for signal-procesing filters -** Progammer: Kevin Rosenberg +** Progammer: Kevin Rosenberg ** Date Started: Aug 1984 ** ** This is part of the CTSim program @@ -30,7 +30,7 @@ int SignalFilter::N_INTEGRAL=500; //static member const int SignalFilter::FILTER_INVALID = -1 ; -const int SignalFilter::FILTER_ABS_BANDLIMIT = 0; // filter times |x| +const int SignalFilter::FILTER_ABS_BANDLIMIT = 0; // filter times |x| const int SignalFilter::FILTER_ABS_G_HAMMING = 1; const int SignalFilter::FILTER_ABS_HANNING = 2; const int SignalFilter::FILTER_ABS_COSINE = 3; @@ -100,13 +100,13 @@ const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const * * SYNOPSIS * f = SignalFilter (filt_type, bw, filterMin, filterMax, n, param, domain, analytic) -* double f Generated filter vector -* int filt_type Type of filter wanted -* double bw Bandwidth of filter -* double filterMin, filterMax Filter limits -* int nFilterPoints Number of points in signal -* double param General input parameter to filters -* int domain FREQUENCY or SPATIAL domain wanted +* double f Generated filter vector +* int filt_type Type of filter wanted +* double bw Bandwidth of filter +* double filterMin, filterMax Filter limits +* int nFilterPoints Number of points in signal +* double param General input parameter to filters +* int domain FREQUENCY or SPATIAL domain wanted */ SignalFilter::SignalFilter (const char* szFilterName, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const char* szDomainName) @@ -140,7 +140,7 @@ SignalFilter::SignalFilter (const char* szFilterName, const char* szDomainName, { m_nFilterPoints = 0; m_dBandwidth = dBandwidth; - m_dFilterParam = dFilterParam; + m_dFilterParam = dFilterParam; m_idFilter = convertFilterNameToID (szFilterName); if (m_idFilter == FILTER_INVALID) { m_fail = true; @@ -173,18 +173,18 @@ SignalFilter::init (const int idFilter, double dFilterMinimum, double dFilterMax m_failMessage = " less than 2"; return; } - + m_nameFilter = convertFilterIDToName (m_idFilter); m_nameDomain = convertDomainIDToName (m_idDomain); m_nFilterPoints = nFilterPoints; - m_dFilterParam = dFilterParam; + m_dFilterParam = dFilterParam; m_dBandwidth = dBandwidth; m_dFilterMin = dFilterMinimum; m_dFilterMax = dFilterMaximum; - + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); m_adFilter = new double [m_nFilterPoints]; - + if (m_idDomain == DOMAIN_FREQUENCY) createFrequencyFilter (m_adFilter); else if (m_idDomain == DOMAIN_SPATIAL) @@ -216,7 +216,7 @@ SignalFilter::createSpatialFilter (double* adFilter) const int center = (m_nFilterPoints - 1) / 2; int sidelen = center; m_adFilter[center] = 4. / (a * a); - + for (int i = 1; i <= sidelen; i++ ) m_adFilter [center + i] = m_adFilter [center - i] = c / (4 * (i * i) - 1); } else { @@ -234,13 +234,13 @@ int 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; } - + return (filterID); } @@ -248,10 +248,10 @@ const char * SignalFilter::convertFilterIDToName (const int filterID) { static const char *name = ""; - + if (filterID >= 0 && filterID < s_iFilterCount) return (s_aszFilterName [filterID]); - + return (name); } @@ -259,10 +259,10 @@ const char * SignalFilter::convertFilterIDToTitle (const int filterID) { static const char *title = ""; - + if (filterID >= 0 && filterID < s_iFilterCount) return (s_aszFilterTitle [filterID]); - + return (title); } @@ -270,13 +270,13 @@ int SignalFilter::convertDomainNameToID (const char* const domainName) { int dID = DOMAIN_INVALID; - + for (int i = 0; i < s_iDomainCount; i++) if (strcasecmp (domainName, s_aszDomainName[i]) == 0) { dID = i; break; } - + return (dID); } @@ -284,10 +284,10 @@ const char * SignalFilter::convertDomainIDToName (const int domainID) { static const char *name = ""; - + if (domainID >= 0 && domainID < s_iDomainCount) return (s_aszDomainName [domainID]); - + return (name); } @@ -295,10 +295,10 @@ const char * SignalFilter::convertDomainIDToTitle (const int domainID) { static const char *title = ""; - + if (domainID >= 0 && domainID < s_iDomainCount) return (s_aszDomainTitle [domainID]); - + return (title); } @@ -307,17 +307,17 @@ double SignalFilter::response (double x) { double response = 0; - + if (m_idDomain == DOMAIN_SPATIAL) response = spatialResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam); else if (m_idDomain == DOMAIN_FREQUENCY) response = frequencyResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam); - + return (response); } -double +double SignalFilter::spatialResponse (int filterID, double bw, double x, double param) { if (haveAnalyticSpatial(filterID)) @@ -331,37 +331,37 @@ SignalFilter::copyFilterData (double* pdFilter, const int iStart, const int nPoi { int iFirst = clamp (iStart, 0, m_nFilterPoints - 1); int iLast = clamp (iFirst + nPoints - 1, 0, m_nFilterPoints - 1); - + for (int i = iFirst; i <= iLast; i++) pdFilter[i - iFirst] = m_adFilter[i]; } /* NAME -* filter_spatial_response_calc Calculate filter by discrete inverse fourier -* transform of filters's frequency -* response +* filter_spatial_response_calc Calculate filter by discrete inverse fourier +* transform of filters's frequency +* response * * SYNOPSIS * y = filter_spatial_response_calc (filt_type, x, m_bw, param, n) -* double y Filter's response in spatial domain -* int filt_type Type of filter (definitions in ct.h) -* double x Spatial position to evaluate filter -* double m_bw Bandwidth of window -* double param General parameter for various filters -* int n Number of points to calculate integrations +* double y Filter's response in spatial domain +* int filt_type Type of filter (definitions in ct.h) +* double x Spatial position to evaluate filter +* double m_bw Bandwidth of window +* double param General parameter for various filters +* int n Number of points to calculate integrations */ -double +double SignalFilter::spatialResponseCalc (double x) const { return (spatialResponseCalc (m_idFilter, m_dBandwidth, x, m_dFilterParam, N_INTEGRAL)); } -double +double SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n) { double zmin, zmax; - + if (filterID == FILTER_TRIANGLE) { zmin = 0; zmax = bw; @@ -370,45 +370,45 @@ SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double par zmax = bw / 2; } double zinc = (zmax - zmin) / (n - 1); - + double z = zmin; double* q = new double [n]; for (int i = 0; i < n; i++, z += zinc) q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x); - + double y = 2 * integrateSimpson (zmin, zmax, q, n); delete q; - + return (y); } /* NAME -* filter_frequency_response Return filter frequency response +* filter_frequency_response Return filter frequency response * * SYNOPSIS * h = filter_frequency_response (filt_type, u, m_bw, param) -* double h Filters frequency response at u -* int filt_type Type of filter -* double u Frequency to evaluate filter at -* double m_bw Bandwidth of filter -* double param General input parameter for various filters +* double h Filters frequency response at u +* int filt_type Type of filter +* double u Frequency to evaluate filter at +* double m_bw Bandwidth of filter +* double param General input parameter for various filters */ -double +double SignalFilter::frequencyResponse (double u) const { return frequencyResponse (m_idFilter, m_dBandwidth, u, m_dFilterParam); } -double +double SignalFilter::frequencyResponse (int filterID, double bw, double u, double param) { double q; double au = fabs (u); double abw = fabs (bw); - + switch (filterID) { case FILTER_BANDLIMIT: if (au >= (abw / 2) + F_EPSILON) @@ -450,7 +450,7 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param q = au * abw * sinc (PI * abw * au, 1.); break; case FILTER_HANNING: - param = 0.5; + param = 0.5; // follow through to G_HAMMING case FILTER_G_HAMMING: if (au >= (abw / 2) + F_EPSILON) @@ -472,27 +472,27 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID); break; } - + return (q); } /* NAME -* filter_spatial_response_analytic Calculate filter by analytic inverse fourier -* transform of filters's frequency -* response +* filter_spatial_response_analytic Calculate filter by analytic inverse fourier +* transform of filters's frequency +* response * * SYNOPSIS * y = filter_spatial_response_analytic (filt_type, x, m_bw, param) -* double y Filter's response in spatial domain -* int filt_type Type of filter (definitions in ct.h) -* double x Spatial position to evaluate filter -* double m_bw Bandwidth of window -* double param General parameter for various filters +* double y Filter's response in spatial domain +* int filt_type Type of filter (definitions in ct.h) +* double x Spatial position to evaluate filter +* double m_bw Bandwidth of window +* double param General parameter for various filters */ -double +double SignalFilter::spatialResponseAnalytic (double x) const { return spatialResponseAnalytic (m_idFilter, m_dBandwidth, x, m_dFilterParam); @@ -502,7 +502,7 @@ const bool SignalFilter::haveAnalyticSpatial (int filterID) { bool haveAnalytic = false; - + switch (filterID) { case FILTER_BANDLIMIT: case FILTER_TRIANGLE: @@ -520,11 +520,11 @@ SignalFilter::haveAnalyticSpatial (int filterID) default: break; } - + return (haveAnalytic); } -double +double SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param) { double q, temp; @@ -532,7 +532,7 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double double w = bw / 2; double b = PI / bw; double b2 = TWOPI / bw; - + switch (filterID) { case FILTER_BANDLIMIT: q = bw * sinc(u * w, 1.0); @@ -581,7 +581,7 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double q = 0; break; } - + return (q); } @@ -590,16 +590,16 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double // Functions that are inline in filter.h -// sinc Return sin(x)/x function +// sinc Return sin(x)/x function // v = sinc (x, mult) // Calculates sin(x * mult) / x; -// integral_abscos Returns integral of u*cos(u) +// integral_abscos Returns integral of u*cos(u) // // q = integral_abscos (u, w) -// double q Integral value -// double u Integration variable -// double w Upper integration boundary +// double q Integral value +// double u Integration variable +// double w Upper integration boundary // Returns the value of integral of u*cos(u)*dV for V = 0 to w