X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Ffilter.cpp;h=397f8a204120926ed04516653b5e1b2ca1d0f654;hp=a88fad5a51c3cb1a663908a81b22dc95f8f7d07b;hb=e8462f7431582627e44906239077f1c696eefba1;hpb=e36dfad3f0818b4c3457fbe7277faa6f4ca28dfe diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index a88fad5..397f8a2 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -1,15 +1,13 @@ /***************************************************************************** ** 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 -** Copyright (C) 1983-2000 Kevin Rosenberg -** -** $Id: filter.cpp,v 1.25 2000/08/19 22:59:06 kevin Exp $ +** Copyright (c) 1983-2009 Kevin Rosenberg ** ** 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 +28,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 +79,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,17 +128,17 @@ 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; - m_dFilterParam = dFilterParam; + m_dFilterParam = dFilterParam; m_idFilter = convertFilterNameToID (szFilterName); if (m_idFilter == FILTER_INVALID) { m_fail = true; @@ -169,7 +175,7 @@ SignalFilter::init (const int idFilter, double dFilterMinimum, double dFilterMax 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; @@ -178,15 +184,15 @@ SignalFilter::init (const int idFilter, double dFilterMinimum, double dFilterMax 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 +214,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); } } } @@ -227,22 +233,23 @@ SignalFilter::convertFilterNameToID (const char *filterName) { int filterID = FILTER_INVALID; - for (int i = 0; i < s_iFilterCount; i++) + 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,24 +258,24 @@ const char * SignalFilter::convertFilterIDToTitle (const int filterID) { static const char *title = ""; - - if (filterID >= 0 && filterID < s_iFilterCount) - return (s_aszFilterTitle [filterID]); + 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) { + for (int i = 0; i < s_iDomainCount; i++) { + if (strcasecmp (domainName, s_aszDomainName[i]) == 0) { dID = i; break; } - + } return (dID); } @@ -278,7 +285,7 @@ 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); } @@ -289,7 +296,7 @@ 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); } @@ -309,7 +316,7 @@ SignalFilter::response (double x) } -double +double SignalFilter::spatialResponse (int filterID, double bw, double x, double param) { if (haveAnalyticSpatial(filterID)) @@ -321,35 +328,35 @@ 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); + 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]; + 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 - */ - -double +* 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 { 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; @@ -364,115 +371,127 @@ SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double par 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 - */ - -double +* 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 { 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 >= 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 - */ - -double +* 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 { return spatialResponseAnalytic (m_idFilter, m_dBandwidth, x, m_dFilterParam); @@ -488,9 +507,11 @@ SignalFilter::haveAnalyticSpatial (int filterID) 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; @@ -502,7 +523,7 @@ SignalFilter::haveAnalyticSpatial (int filterID) return (haveAnalytic); } -double +double SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double param) { double q, temp; @@ -522,6 +543,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 +555,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)); @@ -553,35 +580,26 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double q = 0; break; } - + return (q); } -/* 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 - */