X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=59e02d7f2ba76fbdaa598f9dc53538a988918501;hp=a88fad5a51c3cb1a663908a81b22dc95f8f7d07b;hb=87f312d59cabca5080b481a20314601ea476c4be;hpb=e36dfad3f0818b4c3457fbe7277faa6f4ca28dfe diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index a88fad5..59e02d7 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -7,9 +7,9 @@ ** Date Started: Aug 1984 ** ** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg +** Copyright (c) 1983-2000 Kevin Rosenberg ** -** $Id: filter.cpp,v 1.25 2000/08/19 22:59:06 kevin Exp $ +** $Id: filter.cpp,v 1.40 2003/07/04 21:39:40 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,49 @@ 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_SINC = 1; -const int SignalFilter::FILTER_ABS_G_HAMMING = 2; +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; -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_sinc"}, - {"abs_hamming"}, - {"abs_cosine"}, - {"shepp"}, - {"bandlimit"}, - {"sinc"}, - {"hamming"}, - {"cosine"}, - {"triangle"}, +const int SignalFilter::FILTER_ABS_SINC = 4; +const int SignalFilter::FILTER_SHEPP = 5; +const int SignalFilter::FILTER_BANDLIMIT = 6; +const int SignalFilter::FILTER_SINC = 7; +const int SignalFilter::FILTER_G_HAMMING = 8; +const int SignalFilter::FILTER_HANNING = 9; +const int SignalFilter::FILTER_COSINE = 10; +const int SignalFilter::FILTER_TRIANGLE = 11; + +const int SignalFilter::s_iReconstructFilterCount = 4; + +const char* const SignalFilter::s_aszFilterName[] = { + "abs_bandlimit", + "abs_hamming", + "abs_hanning", + "abs_cosine", + "shepp", + "abs_sinc", + "bandlimit", + "sinc", + "hamming", + "hanning", + "cosine", + "triangle" }; -const char* SignalFilter::s_aszFilterTitle[] = { - {"Abs(w) * Bandlimit"}, - {"Abs(w) * Sinc"}, - {"Abs(w) * Hamming"}, - {"Abs(w) * Cosine"}, - {"Shepp"}, - {"Bandlimit"}, - {"Sinc"}, - {"Hamming"}, - {"Cosine"}, - {"Triangle"}, +const char* const SignalFilter::s_aszFilterTitle[] = { + "Abs(w) * Bandlimit", + "Abs(w) * Hamming", + "Abs(w) * Hanning", + "Abs(w) * Cosine", + "Shepp", + "Abs(w) * Sinc", + "Bandlimit", + "Sinc", + "Hamming", + "Hanning", + "Cosine", + "Triangle" }; const int SignalFilter::s_iFilterCount = sizeof(s_aszFilterName) / sizeof(const char*); @@ -73,36 +81,36 @@ const int SignalFilter::s_iFilterCount = sizeof(s_aszFilterName) / sizeof(const 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* const SignalFilter::s_aszDomainName[] = { + "frequency", + "spatial", }; -const char* SignalFilter::s_aszDomainTitle[] = { - {"Frequency"}, - {"Spatial"}, +const char* const SignalFilter::s_aszDomainTitle[] = { + "Frequency", + "Spatial", }; const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const char*); /* NAME - * SignalFilter::SignalFilter Construct a signal - * - * 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 - */ +* SignalFilter::SignalFilter Construct a signal +* +* 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 +*/ 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_adFilter(NULL), m_fail(false) { m_idFilter = convertFilterNameToID (szFilterName); if (m_idFilter == FILTER_INVALID) { @@ -122,13 +130,13 @@ SignalFilter::SignalFilter (const char* szFilterName, double dFilterMinimum, dou } SignalFilter::SignalFilter (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain) - : m_adFilter(NULL), m_fail(false) +: 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_adFilter(NULL), m_fail(false) { m_nFilterPoints = 0; m_dBandwidth = dBandwidth; @@ -165,7 +173,7 @@ 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; @@ -173,20 +181,20 @@ SignalFilter::init (const int idFilter, double dFilterMinimum, double dFilterMax 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); + createFrequencyFilter (m_adFilter); else if (m_idDomain == DOMAIN_SPATIAL) - createSpatialFilter (m_adFilter); + createSpatialFilter (m_adFilter); } SignalFilter::~SignalFilter (void) { - delete [] m_adFilter; + delete [] m_adFilter; } void @@ -208,16 +216,16 @@ 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 { double x = m_dFilterMin; for (int i = 0; i < m_nFilterPoints; i++, x += m_dFilterInc) { if (haveAnalyticSpatial(m_idFilter)) - m_adFilter[i] = spatialResponseAnalytic (x); + m_adFilter[i] = spatialResponseAnalytic (x); else - m_adFilter[i] = spatialResponseCalc (x); + m_adFilter[i] = spatialResponseCalc (x); } } } @@ -226,24 +234,24 @@ 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); + + return (filterID); } const char * SignalFilter::convertFilterIDToName (const int filterID) { static const char *name = ""; - + if (filterID >= 0 && filterID < s_iFilterCount) - return (s_aszFilterName [filterID]); - + return (s_aszFilterName [filterID]); + return (name); } @@ -251,35 +259,35 @@ const char * SignalFilter::convertFilterIDToTitle (const int filterID) { static const char *title = ""; - + if (filterID >= 0 && filterID < s_iFilterCount) - return (s_aszFilterTitle [filterID]); - + 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) { + if (strcasecmp (domainName, s_aszDomainName[i]) == 0) { dID = i; break; } - - return (dID); + + return (dID); } const char * SignalFilter::convertDomainIDToName (const int domainID) { static const char *name = ""; - + if (domainID >= 0 && domainID < s_iDomainCount) - return (s_aszDomainName [domainID]); - + return (s_aszDomainName [domainID]); + return (name); } @@ -287,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 (s_aszDomainTitle [domainID]); + return (title); } @@ -299,12 +307,12 @@ 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); } @@ -321,27 +329,27 @@ SignalFilter::spatialResponse (int filterID, double bw, double x, double param) 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]; + 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 - * - * 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 - */ +* 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 SignalFilter::spatialResponseCalc (double x) const @@ -353,7 +361,7 @@ double SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double param, int n) { double zmin, zmax; - + if (filterID == FILTER_TRIANGLE) { zmin = 0; zmax = bw; @@ -362,29 +370,30 @@ SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double par zmax = bw / 2; } 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] = 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 - * - * 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 - */ +* 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 SignalFilter::frequencyResponse (double u) const @@ -398,79 +407,90 @@ 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 >= bw / 2) + if (au >= (abw / 2) + F_EPSILON) q = 0.; else q = 1; break; case FILTER_ABS_BANDLIMIT: - if (au >= bw / 2) + if (au >= (abw / 2) + F_EPSILON) q = 0.; else q = au; break; case FILTER_TRIANGLE: - if (au >= bw) + if (au >= (abw / 2) + F_EPSILON) q = 0; else - q = 1 - au / bw; + q = 1 - au / abw; break; case FILTER_COSINE: - if (au >= bw / 2) + if (au >= (abw / 2) + F_EPSILON) q = 0; else - q = cos(PI * u / bw); + q = cos(PI * au / abw); break; case FILTER_ABS_COSINE: - if (au >= bw / 2) + if (au >= (abw / 2) + F_EPSILON) q = 0; else - q = au * cos(PI * u / bw); + q = au * cos(PI * au / abw); break; case FILTER_SINC: - q = bw * sinc (PI * bw * u, 1.); + q = abw * sinc (PI * abw * au, 1.); break; case FILTER_ABS_SINC: - q = au * bw * sinc (PI * bw * u, 1.); + if (au >= (abw / 2) + F_EPSILON) + q = 0; + else + q = au * abw * sinc (PI * abw * au, 1.); break; + case FILTER_HANNING: + param = 0.5; + // follow through to G_HAMMING case FILTER_G_HAMMING: - if (au >= bw / 2) + if (au >= (abw / 2) + F_EPSILON) q = 0; else - q = param + (1 - param) * cos (TWOPI * u / bw); + q = param + (1 - param) * cos (TWOPI * au / abw); break; + case FILTER_ABS_HANNING: + param = 0.5; + // follow through to ABS_G_HAMMING case FILTER_ABS_G_HAMMING: - if (au >= bw / 2) + if (au >= (abw / 2) + F_EPSILON) q = 0; else - q = au * (param + (1 - param) * cos(TWOPI * u / bw)); + q = au * (param + (1 - param) * cos(TWOPI * au / abw)); break; default: q = 0; 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 - * - * 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 - */ +* 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 SignalFilter::spatialResponseAnalytic (double x) const @@ -482,15 +502,17 @@ 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_HANNING: case FILTER_ABS_BANDLIMIT: case FILTER_ABS_COSINE: case FILTER_ABS_G_HAMMING: + case FILTER_ABS_HANNING: case FILTER_SHEPP: case FILTER_SINC: haveAnalytic = true; @@ -498,7 +520,7 @@ SignalFilter::haveAnalyticSpatial (int filterID) default: break; } - + return (haveAnalytic); } @@ -510,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); @@ -522,6 +544,9 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double case FILTER_COSINE: q = sinc(b-u,w) + sinc(b+u,w); break; + case FILTER_HANNING: + param = 0.5; + // follow through to G_HAMMING case FILTER_G_HAMMING: q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w)); break; @@ -531,6 +556,9 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double case FILTER_ABS_COSINE: q = integral_abscos(b-u,w) + integral_abscos(b+u,w); break; + case FILTER_ABS_HANNING: + param = 0.5; + // follow through to ABS_G_HAMMING case FILTER_ABS_G_HAMMING: q = 2 * param * integral_abscos(u,w) + (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w)); @@ -558,30 +586,21 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double } -/* NAME - * sinc Return sin(x)/x function - * - * SYNOPSIS - * v = sinc (x, mult) - * double v sinc value - * double x, mult - * - * DESCRIPTION - * v = sin(x * mult) / x; - */ +// Functions that are inline in filter.h + + +// sinc Return sin(x)/x function +// v = sinc (x, mult) +// Calculates sin(x * mult) / x; + +// 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 +// Returns the value of integral of u*cos(u)*dV for V = 0 to w -/* NAME - * integral_abscos Returns integral of u*cos(u) - * - * SYNOPSIS - * q = integral_abscos (u, w) - * double q Integral value - * double u Integration variable - * double w Upper integration boundary - * - * DESCRIPTION - * Returns the value of integral of u*cos(u)*dV for V = 0 to w - */