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: procsignal.cpp,v 1.5 2000/08/31 08:38:58 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 ******************************************************************************/
30 // FilterMethod ID/Names
31 const int ProcessSignal::FILTER_METHOD_INVALID = -1;
32 const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0;
33 const int ProcessSignal::FILTER_METHOD_FOURIER = 1;
34 const int ProcessSignal::FILTER_METHOD_FOURIER_TABLE = 2;
35 const int ProcessSignal::FILTER_METHOD_FFT = 3;
37 const int ProcessSignal::FILTER_METHOD_FFTW = 4;
38 const int ProcessSignal::FILTER_METHOD_RFFTW =5 ;
40 const char* ProcessSignal::s_aszFilterMethodName[] = {
50 const char* ProcessSignal::s_aszFilterMethodTitle[] = {
53 {"Fouier Trigometric Table"},
57 {"Real/Half-Complex FFTW"},
60 const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*);
62 // FilterGeneration ID/Names
63 const int ProcessSignal::FILTER_GENERATION_INVALID = -1;
64 const int ProcessSignal::FILTER_GENERATION_DIRECT = 0;
65 const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1;
66 const char* ProcessSignal::s_aszFilterGenerationName[] = {
70 const char* ProcessSignal::s_aszFilterGenerationTitle[] = {
74 const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*);
77 // CLASS IDENTIFICATION
80 ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength)
81 : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
83 m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
84 if (m_idFilterMethod == FILTER_METHOD_INVALID) {
86 m_failMessage = "Invalid filter method name ";
87 m_failMessage += szFilterMethodName;
90 m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName);
91 if (m_idFilterGeneration == FILTER_GENERATION_INVALID) {
93 m_failMessage = "Invalid frequency filter name ";
94 m_failMessage += szFilterGenerationName;
97 m_idFilter = SignalFilter::convertFilterNameToID (szFilterName);
98 if (m_idFilter == SignalFilter::FILTER_INVALID) {
100 m_failMessage = "Invalid Filter name ";
101 m_failMessage += szFilterName;
104 m_idDomain = SignalFilter::convertDomainNameToID (szDomainName);
105 if (m_idDomain == SignalFilter::DOMAIN_INVALID) {
107 m_failMessage = "Invalid domain name ";
108 m_failMessage += szDomainName;
112 init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength);
117 ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength)
119 m_idFilter = idFilter;
120 m_idDomain = idDomain;
121 m_idFilterMethod = idFilterMethod;
122 m_idFilterGeneration = idFilterGeneration;
123 m_idGeometry = iGeometry;
124 m_dFocalLength = dFocalLength;
126 if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
130 m_traceLevel = iTraceLevel;
131 m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
132 m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
133 m_dBandwidth = dBandwidth;
134 m_nSignalPoints = nSignalPoints;
135 m_dSignalInc = dSignalIncrement;
136 m_dFilterParam = dFilterParam;
137 m_iZeropad = iZeropad;
138 m_iPreinterpolationFactor = iPreinterpolationFactor;
140 // scale signalInc/BW to signalInc/2 to adjust for imaginary detector
141 // through origin of phantom, see Kak-Slaney Fig 3.22, for Collinear
142 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
147 if (m_idFilterMethod == FILTER_METHOD_FFT) {
149 m_idFilterMethod = FILTER_METHOD_RFFTW;
152 m_failMessage = "FFT not yet implemented";
157 bool m_bFrequencyFiltering = true;
158 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION)
159 m_bFrequencyFiltering = false;
161 // Spatial-based filtering
162 if (! m_bFrequencyFiltering) {
164 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
165 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
166 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
167 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
168 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
169 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
170 m_adFilter = new double[ m_nFilterPoints ];
171 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
172 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
173 m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1;
174 m_dFilterMin = -1. / (2 * m_dSignalInc);
175 m_dFilterMax = 1. / (2 * m_dSignalInc);
176 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
177 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
178 m_adFilter = new double[ m_nFilterPoints ];
179 double adFrequencyFilter [m_nFilterPoints];
180 filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
181 if (m_traceLevel >= Trace::TRACE_PLOT) {
182 SGPDriver sgpDriver ("Frequency Filter: Natural Order");
186 ezplot.ezset ("title Filter Response: Natural Order");
187 ezplot.addCurve (adFrequencyFilter, m_nFilterPoints);
189 cio_put_str ("Press any key to continue");
193 shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints);
194 if (m_traceLevel >= Trace::TRACE_PLOT) {
195 SGPDriver sgpDriver ("Frequency Filter: Fourier Order");
199 ezplot.ezset ("title Filter Response: Fourier Order");
200 ezplot.addCurve (adFrequencyFilter, m_nFilterPoints);
202 cio_put_str ("Press any key to continue");
205 ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
206 if (m_traceLevel >= Trace::TRACE_PLOT) {
207 SGPDriver sgpDriver ("Inverse Fourier Frequency: Fourier Order");
211 ezplot.ezset ("title Inverse Fourier Frequency: Fourier Order");
212 ezplot.addCurve (m_adFilter, m_nFilterPoints);
214 cio_put_str ("Press any key to continue");
217 shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints);
218 if (m_traceLevel >= Trace::TRACE_PLOT) {
219 SGPDriver sgpDriver ("Inverse Fourier Frequency: Natural Order");
223 ezplot.ezset ("title Inverse Fourier Frequency: Natural Order");
224 ezplot.addCurve (m_adFilter, m_nFilterPoints);
226 cio_put_str ("Press any key to continue");
229 for (int i = 0; i < m_nFilterPoints; i++) {
230 m_adFilter[i] /= m_dSignalInc;
233 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
234 for (int i = 0; i < m_nFilterPoints; i++)
235 m_adFilter[i] *= 0.5;
236 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
237 for (int i = 0; i < m_nFilterPoints; i++) {
238 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
239 double sinScale = sin (iDetFromZero * m_dSignalInc);
240 if (fabs(sinScale) < 1E-7)
243 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
244 double dScale = 0.5 * sinScale * sinScale;
245 m_adFilter[i] *= dScale;
248 } // if (spatial filtering)
250 else if (m_bFrequencyFiltering) { // Frequency-based filtering
252 if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
253 // calculate number of filter points with zeropadding
254 m_nFilterPoints = m_nSignalPoints;
255 if (m_iZeropad > 0) {
256 double logBase2 = log(m_nFilterPoints) / log(2);
257 int nextPowerOf2 = static_cast<int>(floor(logBase2));
258 if (logBase2 != floor(logBase2))
260 nextPowerOf2 += (m_iZeropad - 1);
261 m_nFilterPoints = 1 << nextPowerOf2;
262 if (m_traceLevel >= Trace::TRACE_CONSOLE)
263 cout << "nFilterPoints = " << m_nFilterPoints << endl;
265 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
267 if (m_nFilterPoints % 2) { // Odd
268 m_dFilterMin = -1. / (2 * m_dSignalInc);
269 m_dFilterMax = 1. / (2 * m_dSignalInc);
270 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
272 m_dFilterMin = -1. / (2 * m_dSignalInc);
273 m_dFilterMax = 1. / (2 * m_dSignalInc);
274 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints;
275 m_dFilterMax -= m_dFilterInc;
278 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
279 m_adFilter = new double [m_nFilterPoints];
280 filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
282 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
283 for (int i = 0; i < m_nFilterPoints; i++)
284 m_adFilter[i] *= 0.5;
285 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
286 for (int i = 0; i < m_nFilterPoints; i++) {
287 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
288 double sinScale = sin (iDetFromZero * m_dSignalInc);
289 if (fabs(sinScale) < 1E-7)
292 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
293 double dScale = 0.5 * sinScale * sinScale;
294 // m_adFilter[i] *= dScale;
297 if (m_traceLevel >= Trace::TRACE_PLOT) {
298 SGPDriver sgpDriver ("Frequency Filter: Natural Order");
302 ezplot.ezset ("title Filter Filter: Natural Order");
303 ezplot.addCurve (m_adFilter, m_nFilterPoints);
305 cio_put_str ("Press any key to continue");
308 shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
309 if (m_traceLevel >= Trace::TRACE_PLOT) {
310 SGPDriver sgpDriver ("Frequency Filter: Fourier Order");
314 ezplot.ezset ("title Filter Filter: Fourier Order");
315 ezplot.addCurve (m_adFilter, m_nFilterPoints);
317 cio_put_str ("Press any key to continue");
320 } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
321 // calculate number of filter points with zeropadding
322 int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1;
323 m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
324 m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
325 m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1);
326 m_nFilterPoints = nSpatialPoints;
327 if (m_iZeropad > 0) {
328 double logBase2 = log(nSpatialPoints) / log(2);
329 int nextPowerOf2 = static_cast<int>(floor(logBase2));
330 if (logBase2 != floor(logBase2))
332 nextPowerOf2 += (m_iZeropad - 1);
333 m_nFilterPoints = 1 << nextPowerOf2;
335 m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
336 if (m_traceLevel >= Trace::TRACE_CONSOLE)
337 cout << "nFilterPoints = " << m_nFilterPoints << endl;
338 double adSpatialFilter [m_nFilterPoints];
339 SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
340 filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
341 if (m_traceLevel >= Trace::TRACE_PLOT) {
342 SGPDriver sgpDriver ("Spatial Filter: Natural Order");
346 ezplot.ezset ("title Spatial Filter: Natural Order");
347 ezplot.addCurve (adSpatialFilter, nSpatialPoints);
349 cio_put_str ("Press any key to continue");
352 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
353 for (int i = 0; i < m_nFilterPoints; i++)
354 adSpatialFilter[i] *= 0.5;
355 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
356 for (int i = 0; i < m_nFilterPoints; i++) {
357 int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
358 double sinScale = sin (iDetFromZero * m_dSignalInc);
359 if (fabs(sinScale) < 1E-7)
362 sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
363 double dScale = 0.5 * sinScale * sinScale;
364 adSpatialFilter[i] *= dScale;
367 for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
368 adSpatialFilter[i] = 0;
370 m_adFilter = new double [m_nFilterPoints];
371 complex<double> acInverseFilter [m_nFilterPoints];
372 finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
373 for (int i = 0; i < m_nFilterPoints; i++)
374 m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
375 if (m_traceLevel >= Trace::TRACE_PLOT) {
376 SGPDriver sgpDriver ("Spatial Filter: Inverse");
380 ezplot.ezset ("title Spatial Filter: Inverse");
381 ezplot.addCurve (m_adFilter, m_nFilterPoints);
383 cio_put_str ("Press any key to continue");
389 // precalculate sin and cosine tables for fourier transform
390 if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
391 int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1;
392 double angleIncrement = (2. * PI) / m_nFilterPoints;
393 m_adFourierCosTable = new double[ nFourier ];
394 m_adFourierSinTable = new double[ nFourier ];
396 for (int i = 0; i < nFourier; i++) {
397 m_adFourierCosTable[i] = cos (angle);
398 m_adFourierSinTable[i] = sin (angle);
399 angle += angleIncrement;
404 if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) {
405 for (int i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft
406 m_adFilter[i] /= m_nFilterPoints;
409 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
410 m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
411 m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
412 m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
413 m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
414 for (int i = 0; i < m_nFilterPoints; i++)
415 m_adRealFftInput[i] = 0;
416 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
417 m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
418 m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
419 m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
420 m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
421 for (int i = 0; i < m_nFilterPoints; i++)
422 m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
423 for (int i = 0; i < m_nOutputPoints; i++)
424 m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
430 ProcessSignal::~ProcessSignal (void)
432 delete [] m_adFourierSinTable;
433 delete [] m_adFourierCosTable;
434 delete [] m_adFilter;
437 if (m_idFilterMethod == FILTER_METHOD_FFTW) {
438 fftw_destroy_plan(m_complexPlanForward);
439 fftw_destroy_plan(m_complexPlanBackward);
440 delete [] m_adComplexFftInput;
441 delete [] m_adComplexFftSignal;
443 if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
444 rfftw_destroy_plan(m_realPlanForward);
445 rfftw_destroy_plan(m_realPlanBackward);
446 delete [] m_adRealFftInput;
447 delete [] m_adRealFftSignal;
453 ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
455 int fmID = FILTER_METHOD_INVALID;
457 for (int i = 0; i < s_iFilterMethodCount; i++)
458 if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
467 ProcessSignal::convertFilterMethodIDToName (const int fmID)
469 static const char *name = "";
471 if (fmID >= 0 && fmID < s_iFilterMethodCount)
472 return (s_aszFilterMethodName [fmID]);
478 ProcessSignal::convertFilterMethodIDToTitle (const int fmID)
480 static const char *title = "";
482 if (fmID >= 0 && fmID < s_iFilterMethodCount)
483 return (s_aszFilterMethodTitle [fmID]);
490 ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
492 int fgID = FILTER_GENERATION_INVALID;
494 for (int i = 0; i < s_iFilterGenerationCount; i++)
495 if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
504 ProcessSignal::convertFilterGenerationIDToName (const int fgID)
506 static const char *name = "";
508 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
509 return (s_aszFilterGenerationName [fgID]);
515 ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
517 static const char *name = "";
519 if (fgID >= 0 && fgID < s_iFilterGenerationCount)
520 return (s_aszFilterGenerationTitle [fgID]);
526 ProcessSignal::filterSignal (const float constInput[], double output[]) const
528 double input [m_nSignalPoints];
529 for (int i = 0; i < m_nSignalPoints; i++)
530 input[i] = constInput[i];
532 if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
533 for (int i = 0; i < m_nSignalPoints; i++) {
534 int iDetFromCenter = i - (m_nSignalPoints / 2);
535 input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
537 } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
538 for (int i = 0; i < m_nSignalPoints; i++) {
539 int iDetFromCenter = i - (m_nSignalPoints / 2);
540 input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
543 if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
544 for (int i = 0; i < m_nSignalPoints; i++)
545 output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
546 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
547 double inputSignal[m_nFilterPoints];
548 for (int i = 0; i < m_nSignalPoints; i++)
549 inputSignal[i] = input[i];
550 for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
551 inputSignal[i] = 0; // zeropad
552 complex<double> fftSignal[m_nFilterPoints];
553 finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
554 for (int i = 0; i < m_nFilterPoints; i++)
555 fftSignal[i] *= m_adFilter[i];
556 double inverseFourier[m_nFilterPoints];
557 finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);
558 for (int i = 0; i < m_nSignalPoints; i++)
559 output[i] = inverseFourier[i];
560 } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) {
561 double inputSignal[m_nFilterPoints];
562 for (int i = 0; i < m_nSignalPoints; i++)
563 inputSignal[i] = input[i];
564 for (int i = m_nSignalPoints; i < m_nFilterPoints; i++)
565 inputSignal[i] = 0; // zeropad
566 complex<double> fftSignal[m_nFilterPoints];
567 finiteFourierTransform (inputSignal, fftSignal, -1);
568 for (int i = 0; i < m_nFilterPoints; i++)
569 fftSignal[i] *= m_adFilter[i];
570 double inverseFourier[m_nFilterPoints];
571 finiteFourierTransform (fftSignal, inverseFourier, 1);
572 for (int i = 0; i < m_nSignalPoints; i++)
573 output[i] = inverseFourier[i];
576 else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
577 for (int i = 0; i < m_nSignalPoints; i++)
578 m_adRealFftInput[i] = input[i];
580 fftw_real fftOutput [ m_nFilterPoints ];
581 rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
582 for (int i = 0; i < m_nFilterPoints; i++)
583 m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
584 for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
585 m_adRealFftSignal[i] = 0;
587 fftw_real ifftOutput [ m_nOutputPoints ];
588 rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
589 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
590 output[i] = ifftOutput[i];
591 } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
592 for (int i = 0; i < m_nSignalPoints; i++)
593 m_adComplexFftInput[i].re = input[i];
595 fftw_complex fftOutput [ m_nFilterPoints ];
596 fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
597 for (int i = 0; i < m_nFilterPoints; i++) {
598 m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
599 m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
601 fftw_complex ifftOutput [ m_nOutputPoints ];
602 fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
603 for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
604 output[i] = ifftOutput[i].re;
611 * convolve Discrete convolution of two functions
614 * r = convolve (f1, f2, dx, n, np, func_type)
615 * double r Convolved result
616 * double f1[], f2[] Functions to be convolved
617 * double dx Difference between successive x values
618 * int n Array index to center convolution about
619 * int np Number of points in f1 array
620 * int func_type EVEN or ODD or EVEN_AND_ODD function f2
623 * f1 is the projection data, its indices range from 0 to np - 1.
624 * The index for f2, the filter, ranges from -(np-1) to (np-1).
625 * There are 3 ways to handle the negative vertices of f2:
626 * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
627 * All filters used in reconstruction are even.
628 * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
629 * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
630 * for negative indices. Since f2 must range from -(np-1) to (np-1),
631 * if we add (np - 1) to f2's array index, then f2's index will
632 * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
633 * stored at f2[np-1].
637 ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const
641 #if UNOPTIMIZED_CONVOLUTION
642 for (int i = 0; i < np; i++)
643 sum += func[i] * m_adFilter[n - i + (np - 1)];
645 double* f2 = m_adFilter + n + (np - 1);
646 for (int i = 0; i < np; i++)
647 sum += *func++ * *f2--;
655 ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const
659 #if UNOPTIMIZED_CONVOLUTION
660 for (int i = 0; i < np; i++)
661 sum += func[i] * m_adFilter[n - i + (np - 1)];
663 double* f2 = m_adFilter + n + (np - 1);
664 for (int i = 0; i < np; i++)
665 sum += *func++ * *f2--;
673 ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
675 complex<double> complexOutput[n];
677 finiteFourierTransform (input, complexOutput, n, direction);
678 for (int i = 0; i < n; i++)
679 output[i] = complexOutput[i].real();
683 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], const int n, int direction)
690 double angleIncrement = direction * 2 * PI / n;
691 for (int i = 0; i < n; i++) {
694 for (int j = 0; j < n; j++) {
695 double angle = i * j * angleIncrement;
696 sumReal += input[j] * cos(angle);
697 sumImag += input[j] * sin(angle);
703 output[i] = complex<double> (sumReal, sumImag);
709 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
716 double angleIncrement = direction * 2 * PI / n;
717 for (int i = 0; i < n; i++) {
718 complex<double> sum (0,0);
719 for (int j = 0; j < n; j++) {
720 double angle = i * j * angleIncrement;
721 complex<double> exponentTerm (cos(angle), sin(angle));
722 sum += input[j] * exponentTerm;
732 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
739 double angleIncrement = direction * 2 * PI / n;
740 for (int i = 0; i < n; i++) {
742 for (int j = 0; j < n; j++) {
743 double angle = i * j * angleIncrement;
744 sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
753 // Table-based routines
756 ProcessSignal::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
763 for (int i = 0; i < m_nFilterPoints; i++) {
764 double sumReal = 0, sumImag = 0;
765 for (int j = 0; j < m_nFilterPoints; j++) {
766 int tableIndex = i * j;
768 sumReal += input[j] * m_adFourierCosTable[tableIndex];
769 sumImag += input[j] * m_adFourierSinTable[tableIndex];
771 sumReal += input[j] * m_adFourierCosTable[tableIndex];
772 sumImag -= input[j] * m_adFourierSinTable[tableIndex];
776 sumReal /= m_nFilterPoints;
777 sumImag /= m_nFilterPoints;
779 output[i] = complex<double> (sumReal, sumImag);
783 // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
785 ProcessSignal::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
792 for (int i = 0; i < m_nFilterPoints; i++) {
793 double sumReal = 0, sumImag = 0;
794 for (int j = 0; j < m_nFilterPoints; j++) {
795 int tableIndex = i * j;
797 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
798 - input[j].imag() * m_adFourierSinTable[tableIndex];
799 sumImag += input[j].real() * m_adFourierSinTable[tableIndex]
800 + input[j].imag() * m_adFourierCosTable[tableIndex];
802 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
803 - input[j].imag() * -m_adFourierSinTable[tableIndex];
804 sumImag += input[j].real() * -m_adFourierSinTable[tableIndex]
805 + input[j].imag() * m_adFourierCosTable[tableIndex];
809 sumReal /= m_nFilterPoints;
810 sumImag /= m_nFilterPoints;
812 output[i] = complex<double> (sumReal, sumImag);
817 ProcessSignal::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
824 for (int i = 0; i < m_nFilterPoints; i++) {
826 for (int j = 0; j < m_nFilterPoints; j++) {
827 int tableIndex = i * j;
829 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
830 - input[j].imag() * m_adFourierSinTable[tableIndex];
832 sumReal += input[j].real() * m_adFourierCosTable[tableIndex]
833 - input[j].imag() * -m_adFourierSinTable[tableIndex];
837 sumReal /= m_nFilterPoints;
843 // Odd Number of Points
844 // Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
845 // Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
846 // Even Number of Points
847 // Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
848 // Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
851 ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
853 double* pdTemp = new double [n];
855 int iHalfN = (n - 1) / 2;
857 pdTemp[0] = pdVector[iHalfN];
858 for (int i = 0; i < iHalfN; i++)
859 pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
860 for (int i = 0; i < iHalfN; i++)
861 pdTemp[i + iHalfN + 1] = pdVector[i];
864 pdTemp[0] = pdVector[iHalfN];
865 for (int i = 0; i < iHalfN; i++)
866 pdTemp[i + 1] = pdVector[i + iHalfN];
867 for (int i = 0; i < iHalfN - 1; i++)
868 pdTemp[i + iHalfN + 1] = pdVector[i];
871 for (int i = 0; i < n; i++)
872 pdVector[i] = pdTemp[i];
878 ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
880 double* pdTemp = new double [n];
882 int iHalfN = (n - 1) / 2;
884 pdTemp[iHalfN] = pdVector[0];
885 for (int i = 0; i < iHalfN; i++)
886 pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
887 for (int i = 0; i < iHalfN; i++)
888 pdTemp[i] = pdVector[i + iHalfN + 1];
891 pdTemp[iHalfN] = pdVector[0];
892 for (int i = 0; i < iHalfN; i++)
893 pdTemp[i] = pdVector[i + iHalfN];
894 for (int i = 0; i < iHalfN - 1; i++)
895 pdTemp[i + iHalfN + 1] = pdVector[i+1];
898 for (int i = 0; i < n; i++)
899 pdVector[i] = pdTemp[i];