X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=60f97c8e9b7f27319a3fc1a0dcb088f06e360b00;hp=30d47185e3cd08120a9c66c826d6946afb105f71;hb=c00c639073653fac7463a88f2b000f263236550d;hpb=99dd1d6ed10db1f669a5fe6af71225a50fc0ddfb diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index 30d4718..60f97c8 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -1,5 +1,5 @@ /***************************************************************************** -** FILE IDENTIFICATION +** File IDENTIFICATION ** ** Name: filter.cpp ** Purpose: Routines for signal-procesing filters @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: filter.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $ +** $Id: filter.cpp,v 1.33 2001/01/02 16:02:13 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 @@ -27,61 +27,306 @@ #include "ct.h" +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_G_HAMMING = 1; +const int SignalFilter::FILTER_ABS_COSINE = 2; +const int SignalFilter::FILTER_ABS_SINC = 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[] = { + {"abs_bandlimit"}, + {"abs_hamming"}, + {"abs_cosine"}, + {"abs_sinc"}, + {"shepp"}, + {"bandlimit"}, + {"sinc"}, + {"hamming"}, + {"cosine"}, + {"triangle"}, +}; + +const char* SignalFilter::s_aszFilterTitle[] = { + {"Abs(w) * Bandlimit"}, + {"Abs(w) * Hamming"}, + {"Abs(w) * Cosine"}, + {"Abs(w) * Sinc"}, + {"Shepp"}, + {"Bandlimit"}, + {"Sinc"}, + {"Hamming"}, + {"Cosine"}, + {"Triangle"}, +}; + +const int SignalFilter::s_iFilterCount = sizeof(s_aszFilterName) / sizeof(const char*); + + +const int SignalFilter::DOMAIN_INVALID = -1; +const int SignalFilter::DOMAIN_FREQUENCY = 0; +const int SignalFilter::DOMAIN_SPATIAL = 1; + +const char* SignalFilter::s_aszDomainName[] = { + {"frequency"}, + {"spatial"}, +}; + +const char* SignalFilter::s_aszDomainTitle[] = { + {"Frequency"}, + {"Spatial"}, +}; + +const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const char*); + /* 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, 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 * -filter_generate (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint) +SignalFilter::SignalFilter (const char* szFilterName, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const char* szDomainName) + : m_adFilter(NULL), m_fail(false) +{ + m_idFilter = convertFilterNameToID (szFilterName); + if (m_idFilter == FILTER_INVALID) { + m_fail = true; + m_failMessage = "Invalid Filter name "; + m_failMessage += szFilterName; + return; + } + m_idDomain = convertDomainNameToID (szDomainName); + if (m_idDomain == DOMAIN_INVALID) { + m_fail = true; + m_failMessage = "Invalid domain name "; + m_failMessage += szDomainName; + return; + } + init (m_idFilter, dFilterMinimum, dFilterMaximum, nFilterPoints, dBandwidth, dFilterParam, m_idDomain); +} + +SignalFilter::SignalFilter (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain) + : m_adFilter(NULL), m_fail(false) +{ + init (idFilter, dFilterMinimum, dFilterMaximum, nFilterPoints, dBandwidth, dFilterParam, idDomain); +} + +SignalFilter::SignalFilter (const char* szFilterName, const char* szDomainName, double dBandwidth, double dFilterParam) + : m_adFilter(NULL), m_fail(false) +{ + m_nFilterPoints = 0; + m_dBandwidth = dBandwidth; + m_dFilterParam = dFilterParam; + m_idFilter = convertFilterNameToID (szFilterName); + if (m_idFilter == FILTER_INVALID) { + m_fail = true; + m_failMessage = "Invalid Filter name "; + m_failMessage += szFilterName; + return; + } + m_idDomain = convertDomainNameToID (szDomainName); + if (m_idDomain == DOMAIN_INVALID) { + m_fail = true; + m_failMessage = "Invalid domain name "; + m_failMessage += szDomainName; + return; + } +} + +void +SignalFilter::init (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain) { - double *f = new double [n]; - double xinc = (xmax - xmin) / (n - 1); + m_idFilter = idFilter; + m_idDomain = idDomain; + if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) { + m_fail = true; + return; + } + if (nFilterPoints < 2) { + m_fail = true; + m_failMessage = "Number of filter points "; + m_failMessage += nFilterPoints; + m_failMessage = " less than 2"; + return; + } + + m_nameFilter = convertFilterIDToName (m_idFilter); + m_nameDomain = convertDomainIDToName (m_idDomain); + m_nFilterPoints = nFilterPoints; + 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) + createSpatialFilter (m_adFilter); +} + - if (filt_type == FILTER_SHEPP) { - double a = 2 * bw; +SignalFilter::~SignalFilter (void) +{ + delete [] m_adFilter; +} + +void +SignalFilter::createFrequencyFilter (double* adFilter) const +{ + double x; + int i; + for (x = m_dFilterMin, i = 0; i < m_nFilterPoints; x += m_dFilterInc, i++) + adFilter[i] = frequencyResponse (x); +} + + +void +SignalFilter::createSpatialFilter (double* adFilter) const +{ + if (m_idFilter == FILTER_SHEPP) { + double a = 2 * m_dBandwidth; double c = - 4. / (a * a); - int center = (n - 1) / 2; + int center = (m_nFilterPoints - 1) / 2; int sidelen = center; - f[center] = 4. / (a * a); - + m_adFilter[center] = 4. / (a * a); + for (int i = 1; i <= sidelen; i++ ) - f [center + i] = f [center - i] = c / (4 * (i * i) - 1); - } else if (domain == D_FREQ) { - double x; - int i; - for (x = xmin, i = 0; i < n; x += xinc, i++) - f[i] = filter_frequency_response (filt_type, x, bw, param); - } else if (domain == D_SPATIAL) { - double x; - int i; - for (x = xmin, i = 0; i < n; x += xinc, i++) - if (numint == 0) - f[i] = filter_spatial_response_analytic (filt_type, x, bw, param); - else - f[i] = filter_spatial_response_calc (filt_type, x, bw, param, numint); + m_adFilter [center + i] = m_adFilter [center - i] = c / (4 * (i * i) - 1); } else { - sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain); - return (NULL); + double x = m_dFilterMin; + for (int i = 0; i < m_nFilterPoints; i++, x += m_dFilterInc) { + if (haveAnalyticSpatial(m_idFilter)) + m_adFilter[i] = spatialResponseAnalytic (x); + else + m_adFilter[i] = spatialResponseCalc (x); + } } - - return (f); } +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); +} + +const char * +SignalFilter::convertFilterIDToName (const int filterID) +{ + static const char *name = ""; + + if (filterID >= 0 && filterID < s_iFilterCount) + return (s_aszFilterName [filterID]); + + return (name); +} + +const char * +SignalFilter::convertFilterIDToTitle (const int filterID) +{ + static const char *title = ""; + + if (filterID >= 0 && filterID < s_iFilterCount) + return (s_aszFilterTitle [filterID]); + + return (title); +} + +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); +} + +const char * +SignalFilter::convertDomainIDToName (const int domainID) +{ + static const char *name = ""; + + if (domainID >= 0 && domainID < s_iDomainCount) + return (s_aszDomainName [domainID]); + + return (name); +} + +const char * +SignalFilter::convertDomainIDToTitle (const int domainID) +{ + static const char *title = ""; + + if (domainID >= 0 && domainID < s_iDomainCount) + return (s_aszDomainTitle [domainID]); + + return (title); +} + + +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 +SignalFilter::spatialResponse (int filterID, double bw, double x, double param) +{ + if (haveAnalyticSpatial(filterID)) + return spatialResponseAnalytic (filterID, bw, x, param); + else + return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL); +} + +void +SignalFilter::copyFilterData (double* pdFilter, const int iStart, const int nPoints) const +{ + 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 @@ -89,21 +334,27 @@ filter_generate (const FilterType filt_type, double bw, double xmin, double xmax * response * * SYNOPSIS - * y = filter_spatial_response_calc (filt_type, x, bw, param, n) + * 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 bw Bandwidth of window + * double m_bw Bandwidth of window * double param General parameter for various filters * int n Number of points to calculate integrations */ double -filter_spatial_response_calc (int filt_type, double x, double bw, double param, int n) +SignalFilter::spatialResponseCalc (double x) const +{ + return (spatialResponseCalc (m_idFilter, m_dBandwidth, x, m_dFilterParam, N_INTEGRAL)); +} + +double +SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n) { double zmin, zmax; - if (filt_type == FILTER_TRIANGLE) { + if (filterID == FILTER_TRIANGLE) { zmin = 0; zmax = bw; } else { @@ -113,12 +364,13 @@ filter_spatial_response_calc (int filt_type, double x, double bw, double param, double zinc = (zmax - zmin) / (n - 1); double z = zmin; - double q [n]; + double* q = new double [n]; for (int i = 0; i < n; i++, z += zinc) - q[i] = filter_frequency_response (filt_type, z, bw, param) * cos (TWOPI * z * x); + q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x); double y = 2 * integrateSimpson (zmin, zmax, q, n); - + delete q; + return (y); } @@ -127,47 +379,54 @@ filter_spatial_response_calc (int filt_type, double x, double bw, double param, * filter_frequency_response Return filter frequency response * * SYNOPSIS - * h = filter_frequency_response (filt_type, u, bw, param) + * 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 bw Bandwidth of filter + * double m_bw Bandwidth of filter * double param General input parameter for various filters */ double -filter_frequency_response (int filt_type, double u, double bw, double param) +SignalFilter::frequencyResponse (double u) const +{ + return frequencyResponse (m_idFilter, m_dBandwidth, u, m_dFilterParam); +} + + +double +SignalFilter::frequencyResponse (int filterID, double bw, double u, double param) { double q; double au = fabs (u); - switch (filt_type) { + switch (filterID) { case FILTER_BANDLIMIT: - if (au >= bw / 2) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0.; else q = 1; break; case FILTER_ABS_BANDLIMIT: - if (au >= bw / 2) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0.; else q = au; break; case FILTER_TRIANGLE: - if (au >= bw) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = 1 - au / bw; break; case FILTER_COSINE: - if (au >= bw / 2) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = cos(PI * u / bw); break; case FILTER_ABS_COSINE: - if (au >= bw / 2) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = au * cos(PI * u / bw); @@ -179,22 +438,20 @@ filter_frequency_response (int filt_type, double u, double bw, double param) q = au * bw * sinc (PI * bw * u, 1.); break; case FILTER_G_HAMMING: - if (au >= bw / 2) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = param + (1 - param) * cos (TWOPI * u / bw); break; case FILTER_ABS_G_HAMMING: - if (au >= bw / 2) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = au * (param + (1 - param) * cos(TWOPI * u / bw)); break; default: q = 0; - sys_error (ERR_WARNING, - "Frequency response for filter %d not implemented [filter_frequency_response]", - filt_type); + sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID); break; } return (q); @@ -208,16 +465,46 @@ filter_frequency_response (int filt_type, double u, double bw, double param) * response * * SYNOPSIS - * y = filter_spatial_response_analytic (filt_type, x, bw, param) + * 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 bw Bandwidth of window + * double m_bw Bandwidth of window * double param General parameter for various filters */ double -filter_spatial_response_analytic (int filt_type, double x, double bw, double param) +SignalFilter::spatialResponseAnalytic (double x) const +{ + return spatialResponseAnalytic (m_idFilter, m_dBandwidth, x, m_dFilterParam); +} + +const bool +SignalFilter::haveAnalyticSpatial (int filterID) +{ + bool haveAnalytic = false; + + switch (filterID) { + case FILTER_BANDLIMIT: + case FILTER_TRIANGLE: + case FILTER_COSINE: + case FILTER_G_HAMMING: + case FILTER_ABS_BANDLIMIT: + case FILTER_ABS_COSINE: + case FILTER_ABS_G_HAMMING: + case FILTER_SHEPP: + case FILTER_SINC: + haveAnalytic = true; + break; + default: + break; + } + + return (haveAnalytic); +} + +double +SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param) { double q, temp; double u = TWOPI * x; @@ -225,7 +512,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par double b = PI / bw; double b2 = TWOPI / bw; - switch (filt_type) { + switch (filterID) { case FILTER_BANDLIMIT: q = bw * sinc(u * w, 1.0); break; @@ -237,8 +524,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par q = sinc(b-u,w) + sinc(b+u,w); break; case FILTER_G_HAMMING: - q = 2 * param * sin(u*w)/u + (1-param) * - (sinc(b2-u, w) + sinc(b2+u, w)); + q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w)); break; case FILTER_ABS_BANDLIMIT: q = 2 * integral_abscos (u, w); @@ -264,9 +550,7 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par break; case FILTER_ABS_SINC: default: - sys_error (ERR_WARNING, - "Analytic filter type %d not implemented [filter_spatial_response_analytic]", - filt_type); + sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID); q = 0; break; } @@ -287,12 +571,6 @@ filter_spatial_response_analytic (int filt_type, double x, double bw, double par * v = sin(x * mult) / x; */ -double -sinc (double x, double mult) -{ - return (fabs(x) > F_EPSILON ? (sin (x * mult) / x) : 1.0); -} - /* NAME * integral_abscos Returns integral of u*cos(u) @@ -307,13 +585,4 @@ sinc (double x, double mult) * Returns the value of integral of u*cos(u)*dV for V = 0 to w */ -double -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); -} -