1 /*****************************************************************************
5 ** Purpose: Routines for signal-procesing filters
6 ** Progammer: Kevin Rosenberg
7 ** Date Started: Aug 1984
9 ** This is part of the CTSim program
10 ** Copyright (C) 1983-2000 Kevin Rosenberg
12 ** $Id: filter.cpp,v 1.4 2000/06/25 17:32:24 kevin Exp $
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.
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.
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 ******************************************************************************/
32 * SignalFilter::SignalFilter Construct a signal
35 * f = SignalFilter (filt_type, bw, xmin, xmax, n, param, domain, analytic)
36 * double f Generated filter vector
37 * int filt_type Type of filter wanted
38 * double bw Bandwidth of filter
39 * double xmin, xmax Filter limits
40 * int n Number of points in filter
41 * double param General input parameter to filters
42 * int domain FREQ or SPATIAL domain wanted
43 * int numint Number if intervals for calculating discrete inverse fourier xform
44 * for spatial domain filters. For ANALYTIC solutions, use numint = 0
47 SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numIntegral = 0)
50 m_idFilter = convertFilterNameToID (filterName);
51 if (m_idFilter == FILTER_INVALID) {
53 m_failMessage = "Invalid Filter name ";
54 m_failMessage += filterName;
57 m_idDomain = convertDomainNameToID (domainName);
58 if (m_idDomain == DOMAIN_INVALID) {
60 m_failMessage = "Invalid domain name ";
61 m_failMessage += domainName;
64 init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numIntegral);
67 SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numIntegral = 0)
69 init (filterID, bw, xmin, xmax, n, param, domainID, numIntegral);
72 SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0)
77 m_filterParam = param;
78 m_numIntegral = numIntegral;
79 m_idFilter = convertFilterNameToID (filterName);
80 if (m_idFilter == FILTER_INVALID) {
82 m_failMessage = "Invalid Filter name ";
83 m_failMessage += filterName;
86 m_idDomain = convertDomainNameToID (domainName);
87 if (m_idDomain == DOMAIN_INVALID) {
89 m_failMessage = "Invalid domain name ";
90 m_failMessage += domainName;
96 SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
99 m_idFilter = filterID;
100 m_idDomain = domainID;
101 if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) {
105 m_nameFilter = convertFilterIDToName (m_idFilter);
106 m_nameDomain = convertDomainIDToName (m_idDomain);
111 m_numIntegral = numint;
112 m_filterParam = param;
113 m_vecFilter = new double[n];
115 double xinc = (m_xmax - m_xmin) / (m_nPoints - 1);
117 if (m_idFilter == FILTER_SHEPP) {
119 double c = - 4. / (a * a);
120 int center = (m_nPoints - 1) / 2;
121 int sidelen = center;
122 m_vecFilter[center] = 4. / (a * a);
124 for (int i = 1; i <= sidelen; i++ )
125 m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
126 } else if (m_idDomain == DOMAIN_FREQ) {
129 for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
130 m_vecFilter[i] = frequencyResponse (x, param);
131 } else if (m_idDomain == DOMAIN_SPATIAL) {
134 for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
136 m_vecFilter[i] = spatialResponseAnalytic (x, param);
138 m_vecFilter[i] = spatialResponseCalc (x, param, numint);
140 m_failMessage = "Illegal domain name ";
141 m_failMessage += m_idDomain;
146 SignalFilter::~SignalFilter (void)
152 SignalFilter::FilterID
153 SignalFilter::convertFilterNameToID (const char *filterName)
155 FilterID filterID = FILTER_INVALID;
157 if (strcasecmp (filterName, FILTER_BANDLIMIT_STR) == 0)
158 filterID = FILTER_BANDLIMIT;
159 else if (strcasecmp (filterName, FILTER_HAMMING_STR) == 0)
160 filterID = FILTER_G_HAMMING;
161 else if (strcasecmp (filterName, FILTER_SINC_STR) == 0)
162 filterID = FILTER_SINC;
163 else if (strcasecmp (filterName, FILTER_COS_STR) == 0)
164 filterID = FILTER_COSINE;
165 else if (strcasecmp (filterName, FILTER_TRIANGLE_STR) == 0)
166 filterID = FILTER_TRIANGLE;
167 else if (strcasecmp (filterName, FILTER_ABS_BANDLIMIT_STR) == 0)
168 filterID = FILTER_ABS_BANDLIMIT;
169 else if (strcasecmp (filterName, FILTER_ABS_HAMMING_STR) == 0)
170 filterID = FILTER_ABS_G_HAMMING;
171 else if (strcasecmp (filterName, FILTER_ABS_SINC_STR) == 0)
172 filterID = FILTER_ABS_SINC;
173 else if (strcasecmp (filterName, FILTER_ABS_COS_STR) == 0)
174 filterID = FILTER_ABS_COSINE;
175 else if (strcasecmp (filterName, FILTER_SHEPP_STR) == 0)
176 filterID = FILTER_SHEPP;
182 SignalFilter::convertFilterIDToName (const FilterID filterID)
184 const char *name = "";
186 if (filterID == FILTER_SHEPP)
187 name = FILTER_SHEPP_STR;
188 else if (filterID == FILTER_ABS_COSINE)
189 name = FILTER_ABS_COS_STR;
190 else if (filterID == FILTER_ABS_SINC)
191 name = FILTER_ABS_SINC_STR;
192 else if (filterID == FILTER_ABS_G_HAMMING)
193 name = FILTER_ABS_HAMMING_STR;
194 else if (filterID == FILTER_ABS_BANDLIMIT)
195 name = FILTER_ABS_BANDLIMIT_STR;
196 else if (filterID == FILTER_COSINE)
197 name = FILTER_COS_STR;
198 else if (filterID == FILTER_SINC)
199 name = FILTER_SINC_STR;
200 else if (filterID == FILTER_G_HAMMING)
201 name = FILTER_HAMMING_STR;
202 else if (filterID == FILTER_BANDLIMIT)
203 name = FILTER_BANDLIMIT_STR;
204 else if (filterID == FILTER_TRIANGLE)
205 name = FILTER_TRIANGLE_STR;
210 const SignalFilter::DomainID
211 SignalFilter::convertDomainNameToID (const char* const domainName)
213 DomainID dID = DOMAIN_INVALID;
215 if (strcasecmp (domainName, DOMAIN_SPATIAL_STR) == 0)
216 dID = DOMAIN_SPATIAL;
217 else if (strcasecmp (domainName, DOMAIN_FREQ_STR) == 0)
224 SignalFilter::convertDomainIDToName (const DomainID domain)
226 const char *name = "";
228 if (domain == DOMAIN_SPATIAL)
229 return (DOMAIN_SPATIAL_STR);
230 else if (domain == DOMAIN_FREQ)
231 return (DOMAIN_FREQ_STR);
238 SignalFilter::response (double x)
242 if (m_idDomain == DOMAIN_SPATIAL)
243 response = spatialResponse (m_idFilter, m_bw, x, m_filterParam, m_numIntegral);
244 else if (m_idDomain == DOMAIN_FREQ)
245 response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam);
252 SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double param, int nIntegral = 0)
255 return spatialResponseAnalytic (filterID, bw, x, param);
257 return spatialResponseCalc (filterID, bw, x, param, nIntegral);
261 * filter_spatial_response_calc Calculate filter by discrete inverse fourier
262 * transform of filters's frequency
266 * y = filter_spatial_response_calc (filt_type, x, m_bw, param, n)
267 * double y Filter's response in spatial domain
268 * int filt_type Type of filter (definitions in ct.h)
269 * double x Spatial position to evaluate filter
270 * double m_bw Bandwidth of window
271 * double param General parameter for various filters
272 * int n Number of points to calculate integrations
276 SignalFilter::spatialResponseCalc (double x, double param, int nIntegral) const
278 return (spatialResponseCalc (m_idFilter, m_bw, x, param, nIntegral));
282 SignalFilter::spatialResponseCalc (FilterID filterID, double bw, double x, double param, int n)
286 if (filterID == FILTER_TRIANGLE) {
293 double zinc = (zmax - zmin) / (n - 1);
297 for (int i = 0; i < n; i++, z += zinc)
298 q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
300 double y = 2 * integrateSimpson (zmin, zmax, q, n);
307 * filter_frequency_response Return filter frequency response
310 * h = filter_frequency_response (filt_type, u, m_bw, param)
311 * double h Filters frequency response at u
312 * int filt_type Type of filter
313 * double u Frequency to evaluate filter at
314 * double m_bw Bandwidth of filter
315 * double param General input parameter for various filters
319 SignalFilter::frequencyResponse (double u, double param) const
321 return frequencyResponse (m_idFilter, m_bw, u, param);
326 SignalFilter::frequencyResponse (FilterID filterID, double bw, double u, double param)
329 double au = fabs (u);
332 case FILTER_BANDLIMIT:
338 case FILTER_ABS_BANDLIMIT:
344 case FILTER_TRIANGLE:
354 q = cos(PI * u / bw);
356 case FILTER_ABS_COSINE:
360 q = au * cos(PI * u / bw);
363 q = bw * sinc (PI * bw * u, 1.);
365 case FILTER_ABS_SINC:
366 q = au * bw * sinc (PI * bw * u, 1.);
368 case FILTER_G_HAMMING:
372 q = param + (1 - param) * cos (TWOPI * u / bw);
374 case FILTER_ABS_G_HAMMING:
378 q = au * (param + (1 - param) * cos(TWOPI * u / bw));
382 sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
391 * filter_spatial_response_analytic Calculate filter by analytic inverse fourier
392 * transform of filters's frequency
396 * y = filter_spatial_response_analytic (filt_type, x, m_bw, param)
397 * double y Filter's response in spatial domain
398 * int filt_type Type of filter (definitions in ct.h)
399 * double x Spatial position to evaluate filter
400 * double m_bw Bandwidth of window
401 * double param General parameter for various filters
405 SignalFilter::spatialResponseAnalytic (double x, double param) const
407 return spatialResponseAnalytic (m_idFilter, m_bw, x, param);
411 SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, double param)
414 double u = TWOPI * x;
417 double b2 = TWOPI / bw;
420 case FILTER_BANDLIMIT:
421 q = bw * sinc(u * w, 1.0);
423 case FILTER_TRIANGLE:
424 temp = sinc (u * w, 1.0);
425 q = bw * temp * temp;
428 q = sinc(b-u,w) + sinc(b+u,w);
430 case FILTER_G_HAMMING:
431 q = 2 * param * sin(u*w)/u + (1-param) * (sinc(b2-u, w) + sinc(b2+u, w));
433 case FILTER_ABS_BANDLIMIT:
434 q = 2 * integral_abscos (u, w);
436 case FILTER_ABS_COSINE:
437 q = integral_abscos(b-u,w) + integral_abscos(b+u,w);
439 case FILTER_ABS_G_HAMMING:
440 q = 2 * param * integral_abscos(u,w) +
441 (1-param)*(integral_abscos(u-b2,w)+integral_abscos(u+b2,w));
445 q = 4. / (PI * bw * bw);
447 q = fabs ((2 / bw) * sin (u * w)) * sinc (u * w, 1.) * sinc (u * w, 1.);
450 if (fabs (x) < bw / 2)
455 case FILTER_ABS_SINC:
457 sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
467 * sinc Return sin(x)/x function
471 * double v sinc value
475 * v = sin(x * mult) / x;
480 * integral_abscos Returns integral of u*cos(u)
483 * q = integral_abscos (u, w)
484 * double q Integral value
485 * double u Integration variable
486 * double w Upper integration boundary
489 * Returns the value of integral of u*cos(u)*dV for V = 0 to w
493 SignalFilter::integral_abscos (double u, double w)
495 if (fabs (u) > F_EPSILON)
496 return (cos(u * w) - 1) / (u * u) + w / u * sin (u * w);
503 * convolve Discrete convolution of two functions
506 * r = convolve (f1, f2, dx, n, np, func_type)
507 * double r Convolved result
508 * double f1[], f2[] Functions to be convolved
509 * double dx Difference between successive x values
510 * int n Array index to center convolution about
511 * int np Number of points in f1 array
512 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
515 * f1 is the projection data, its indices range from 0 to np - 1.
516 * The index for f2, the filter, ranges from -(np-1) to (np-1).
517 * There are 3 ways to handle the negative vertices of f2:
518 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
519 * All filters used in reconstruction are even.
520 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
521 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
522 * for negative indices. Since f2 must range from -(np-1) to (np-1),
523 * if we add (np - 1) to f2's array index, then f2's index will
524 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
525 * stored at f2[np-1].
529 SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const
533 #if UNOPTIMIZED_CONVOLUTION
534 for (int i = 0; i < np; i++)
535 sum += func[i] * m_vecFilter[n - i + (np - 1)];
537 double* f2 = m_vecFilter + n + (np - 1);
538 for (int i = 0; i < np; i++)
539 sum += *func++ * *f2--;
547 SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const
551 #if UNOPTIMIZED_CONVOLUTION
552 for (int i = 0; i < np; i++)
553 sum += func[i] * m_vecFilter[n - i + (np - 1)];
555 double* f2 = m_vecFilter + n + (np - 1);
556 for (int i = 0; i < np; i++)
557 sum += *func++ * *f2--;