X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=dd4d40a05dc9b27271aed80e4eeb69d261ddcadb;hp=70d759cbacc471f96ab39b3d8fb3ea3dd990e378;hb=f13a8c004b8f182b42d9e4df2bcd7c7f030bf1ad;hpb=2a39ee3b125e3e2e68bbba2ac15a65039456ff7e diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 70d759c..dd4d40a 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -1,15 +1,13 @@ /***************************************************************************** ** File IDENTIFICATION -** -** Name: filter.cpp -** Purpose: Routines for signal-procesing filters -** Progammer: Kevin Rosenberg -** Date Started: Aug 1984 ** -** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg +** Name: procsignal.cpp +** Purpose: Routines for processing signals and projections +** Progammer: Kevin Rosenberg +** Date Started: Aug 1984 ** -** $Id: procsignal.cpp,v 1.2 2000/08/22 07:02:48 kevin Exp $ +** This is part of the CTSim program +** 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 @@ -27,6 +25,10 @@ #include "ct.h" +#ifdef HAVE_WXWINDOWS +#include "nographics.h" +#endif + // FilterMethod ID/Names const int ProcessSignal::FILTER_METHOD_INVALID = -1; const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0; @@ -37,24 +39,24 @@ const int ProcessSignal::FILTER_METHOD_FFT = 3; const int ProcessSignal::FILTER_METHOD_FFTW = 4; const int ProcessSignal::FILTER_METHOD_RFFTW =5 ; #endif -const char* ProcessSignal::s_aszFilterMethodName[] = { - {"convolution"}, - {"fourier"}, - {"fouier_table"}, - {"fft"}, +const char* const ProcessSignal::s_aszFilterMethodName[] = { + "convolution", + "fourier", + "fouier-table", + "fft", #if HAVE_FFTW - {"fftw"}, - {"rfftw"}, + "fftw", + "rfftw", #endif }; -const char* ProcessSignal::s_aszFilterMethodTitle[] = { - {"Convolution"}, - {"Fourier"}, - {"Fouier Trigometric Table"}, - {"FFT"}, +const char* const ProcessSignal::s_aszFilterMethodTitle[] = { + "Convolution", + "Fourier", + "Fouier Trigometric Table", + "FFT", #if HAVE_FFTW - {"FFTW"}, - {"Real/Half-Complex FFTW"}, + "FFTW", + "Real/Half-Complex FFTW", #endif }; const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*); @@ -63,13 +65,13 @@ const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / const int ProcessSignal::FILTER_GENERATION_INVALID = -1; const int ProcessSignal::FILTER_GENERATION_DIRECT = 0; const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1; -const char* ProcessSignal::s_aszFilterGenerationName[] = { - {"direct"}, - {"inverse_fourier"}, +const char* const ProcessSignal::s_aszFilterGenerationName[] = { + "direct", + "inverse-fourier", }; -const char* ProcessSignal::s_aszFilterGenerationTitle[] = { - {"Direct"}, - {"Inverse Fourier"}, +const char* const ProcessSignal::s_aszFilterGenerationTitle[] = { + "Direct", + "Inverse Fourier", }; const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*); @@ -77,8 +79,11 @@ const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGeneration // CLASS IDENTIFICATION // ProcessSignal // -ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad = 0, int iPreinterpolationFactor = 1, int iTraceLevel = TRACE_NONE) - : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) +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, double dSourceDetectorLength, SGP* pSGP) + : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) { m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName); if (m_idFilterMethod == FILTER_METHOD_INVALID) { @@ -109,17 +114,27 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth return; } - init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel); + init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, + m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, + dSourceDetectorLength, pSGP); } void -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) +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, double dSourceDetectorLength, SGP* pSGP) { + int i; m_idFilter = idFilter; m_idDomain = idDomain; m_idFilterMethod = idFilterMethod; m_idFilterGeneration = idFilterGeneration; + m_idGeometry = iGeometry; + m_dFocalLength = dFocalLength; + m_dSourceDetectorLength = dSourceDetectorLength; + if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) { m_fail = true; return; @@ -130,10 +145,18 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dBandwidth = dBandwidth; m_nSignalPoints = nSignalPoints; m_dSignalInc = dSignalIncrement; - m_dFilterParam = dFilterParam; + m_dFilterParam = dFilterParam; m_iZeropad = iZeropad; m_iPreinterpolationFactor = iPreinterpolationFactor; + // scale signalInc/BW to adjust for imaginary detector through origin of phantom + // see Kak-Slaney Fig 3.22, for Collinear diagram + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { + double dEquilinearScale = m_dSourceDetectorLength / m_dFocalLength; + m_dSignalInc /= dEquilinearScale; + m_dBandwidth *= dEquilinearScale; + } + if (m_idFilterMethod == FILTER_METHOD_FFT) { #if HAVE_FFTW m_idFilterMethod = FILTER_METHOD_RFFTW; @@ -152,131 +175,149 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw if (! m_bFrequencyFiltering) { if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { - m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; - m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); - m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); - m_adFilter = new double[ m_nFilterPoints ]; - filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); + m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; + m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); + m_adFilter = new double[ m_nFilterPoints ]; + filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { - m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); - m_adFilter = new double[ m_nFilterPoints ]; - double adFrequencyFilter [m_nFilterPoints]; - filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints); - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Frequency Filter: Natural Order"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Filter Response: Natural Order"); - ezplot.addCurve (adFrequencyFilter, m_nFilterPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); - } - - shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Frequency Filter: Fourier Order"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Filter Response: Fourier Order"); - ezplot.addCurve (adFrequencyFilter, m_nFilterPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); - } - ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1); - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Inverse Fourier Frequency: Fourier Order"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Inverse Fourier Frequency: Fourier Order"); - ezplot.addCurve (m_adFilter, m_nFilterPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); - } - shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints); - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Inverse Fourier Frequency: Natural Order"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Inverse Fourier Frequency: Natural Order"); - ezplot.addCurve (m_adFilter, m_nFilterPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); - } - for (int i = 0; i < m_nFilterPoints; i++) { - m_adFilter[i] /= m_dSignalInc; - } + m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); + m_adFilter = new double[ m_nFilterPoints ]; + double* adFrequencyFilter = new double [m_nFilterPoints]; + filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints); +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif + Fourier::shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); +#ifdef HAVE_SGP + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Filter Response: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (adFrequencyFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif + ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, FORWARD); + delete adFrequencyFilter; +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif + Fourier::shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints); +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Inverse Fourier Frequency: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif + for (i = 0; i < m_nFilterPoints; i++) { + m_adFilter[i] /= m_dSignalInc; + } } - } + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] *= 0.5; + } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { + for (i = 0; i < m_nFilterPoints; i++) { + int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); + double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc); + double dScale = 0.5 * sinScale * sinScale; + m_adFilter[i] *= dScale; + } +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Scaled Inverse Fourier Frequency: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif + } // if (geometry) + } // if (spatial filtering) - // Frequency-based filtering - else if (m_bFrequencyFiltering) { + else if (m_bFrequencyFiltering) { // Frequency-based filtering if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { // calculate number of filter points with zeropadding - m_nFilterPoints = m_nSignalPoints; - if (m_iZeropad > 0) { - double logBase2 = log(m_nFilterPoints) / log(2); - int nextPowerOf2 = static_cast(floor(logBase2)); - if (logBase2 != floor(logBase2)) - nextPowerOf2++; - nextPowerOf2 += (m_iZeropad - 1); - m_nFilterPoints = 1 << nextPowerOf2; - if (m_traceLevel >= TRACE_TEXT) - cout << "nFilterPoints = " << m_nFilterPoints << endl; - } + m_nFilterPoints = addZeropadFactor (m_nSignalPoints, m_iZeropad); m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - if (m_nFilterPoints % 2) { // Odd - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + if (isOdd (m_nFilterPoints)) { // Odd + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); } else { // Even - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; - m_dFilterMax -= m_dFilterInc; + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; + m_dFilterMax -= m_dFilterInc; } - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, + m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); m_adFilter = new double [m_nFilterPoints]; filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Frequency Filter: Natural Order"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Filter Filter: Natural Order"); - ezplot.addCurve (m_adFilter, m_nFilterPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); + +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Frequency Filter: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); } - shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints); - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Frequency Filter: Fourier Order"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Filter Filter: Fourier Order"); - ezplot.addCurve (m_adFilter, m_nFilterPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); - } +#endif + + // This works fairly well. I'm not sure why since scaling for geometries is done on + // frequency filter rather than spatial filter as it should be. + // It gives values slightly off than freq/inverse filtering + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] *= 0.5; + } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { + for (i = 0; i < m_nFilterPoints; i++) { + int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); + double sinScale = 1 / SignalFilter::sinc (iDetFromZero * m_dSignalInc); + double dScale = 0.5 * sinScale * sinScale; + m_adFilter[i] *= dScale; + } + } +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif + Fourier::shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints); +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Filter Geometry Scaled: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); + } +#endif + + // FILTERING: FREQUENCY - INVERSE FOURIER + } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { // calculate number of filter points with zeropadding int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1; @@ -285,60 +326,83 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1); m_nFilterPoints = nSpatialPoints; if (m_iZeropad > 0) { - double logBase2 = log(nSpatialPoints) / log(2); - int nextPowerOf2 = static_cast(floor(logBase2)); - if (logBase2 != floor(logBase2)) - nextPowerOf2++; - nextPowerOf2 += (m_iZeropad - 1); - m_nFilterPoints = 1 << nextPowerOf2; + double logBase2 = log(nSpatialPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)); + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + nextPowerOf2 += (m_iZeropad - 1); + m_nFilterPoints = 1 << nextPowerOf2; } m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - if (m_traceLevel >= TRACE_TEXT) - cout << "nFilterPoints = " << m_nFilterPoints << endl; - double adSpatialFilter [m_nFilterPoints]; - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); +#if defined(DEBUG) || defined(_DEBUG) + if (m_traceLevel >= Trace::TRACE_CONSOLE) + sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints); +#endif + double* adSpatialFilter = new double [m_nFilterPoints]; + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, + m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints); - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Spatial Filter: Natural Order"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Spatial Filter: Natural Order"); - ezplot.addCurve (adSpatialFilter, nSpatialPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot;; + dlgEZPlot.getEZPlot()->ezset ("title Spatial Filter: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints); + dlgEZPlot.ShowModal(); } - for (int i = nSpatialPoints; i < m_nFilterPoints; i++) - adSpatialFilter[i] = 0; +#endif + + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { + for (i = 0; i < nSpatialPoints; i++) + adSpatialFilter[i] *= 0.5; + } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { + for (i = 0; i < nSpatialPoints; i++) { + int iDetFromZero = i - ((nSpatialPoints - 1) / 2); + double sinScale = sin (iDetFromZero * m_dSignalInc); + if (fabs(sinScale) < 1E-7) + sinScale = 1; + else + sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double dScale = 0.5 * sinScale * sinScale; + adSpatialFilter[i] *= dScale; + } + } +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot;; + dlgEZPlot.getEZPlot()->ezset ("title Scaled Spatial Filter: Natural Order"); + dlgEZPlot.getEZPlot()->addCurve (adSpatialFilter, nSpatialPoints); + dlgEZPlot.ShowModal(); + } +#endif + for (i = nSpatialPoints; i < m_nFilterPoints; i++) + adSpatialFilter[i] = 0; m_adFilter = new double [m_nFilterPoints]; - complex acInverseFilter [m_nFilterPoints]; - finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1); - for (int i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc; - if (m_traceLevel >= TRACE_PLOT) { - SGPDriver sgpDriver ("Spatial Filter: Inverse"); - SGP sgp (sgpDriver); - EZPlot ezplot (sgp); - - ezplot.ezset ("title Spatial Filter: Inverse"); - ezplot.addCurve (m_adFilter, m_nFilterPoints); - ezplot.plot(); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); + std::complex* acInverseFilter = new std::complex [m_nFilterPoints]; + finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, BACKWARD); + delete adSpatialFilter; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] = std::abs (acInverseFilter[i]) * m_dSignalInc; + delete acInverseFilter; +#if defined(HAVE_WXWINDOWS) && (defined(DEBUG) || defined(_DEBUG)) + if (g_bRunningWXWindows && m_traceLevel > 0) { + EZPlotDialog dlgEZPlot; + dlgEZPlot.getEZPlot()->ezset ("title Fourier Scaled Spatial Filter: Fourier Order"); + dlgEZPlot.getEZPlot()->addCurve (m_adFilter, m_nFilterPoints); + dlgEZPlot.ShowModal(); } +#endif } } - + // precalculate sin and cosine tables for fourier transform if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { - int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1; + int nFourier = imax (m_nFilterPoints,m_nOutputPoints) * imax (m_nFilterPoints, m_nOutputPoints) + 1; double angleIncrement = (2. * PI) / m_nFilterPoints; m_adFourierCosTable = new double[ nFourier ]; m_adFourierSinTable = new double[ nFourier ]; double angle = 0; - for (int i = 0; i < nFourier; i++) { + for (i = 0; i < nFourier; i++) { m_adFourierCosTable[i] = cos (angle); m_adFourierSinTable[i] = sin (angle); angle += angleIncrement; @@ -347,50 +411,59 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw #if HAVE_FFTW if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) { - for (int i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft + for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft m_adFilter[i] /= m_nFilterPoints; } if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE); - m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE); - m_adRealFftInput = new fftw_real [ m_nFilterPoints ]; - m_adRealFftSignal = new fftw_real [ m_nOutputPoints ]; - for (int i = 0; i < m_nFilterPoints; i++) + m_adRealFftInput = static_cast(fftw_malloc (sizeof(double) * m_nFilterPoints)); + m_adRealFftOutput = static_cast(fftw_malloc (sizeof(double) * m_nFilterPoints)); + m_realPlanForward = fftw_plan_r2r_1d (m_nFilterPoints, m_adRealFftInput, m_adRealFftOutput, FFTW_R2HC, FFTW_ESTIMATE); + m_adRealFftSignal = static_cast(fftw_malloc (sizeof(double) * m_nOutputPoints)); + m_adRealFftBackwardOutput = static_cast(fftw_malloc (sizeof(double) * m_nOutputPoints)); + m_realPlanBackward = fftw_plan_r2r_1d (m_nOutputPoints, m_adRealFftSignal, m_adRealFftBackwardOutput, FFTW_HC2R, FFTW_ESTIMATE); + for (i = 0; i < m_nFilterPoints; i++) m_adRealFftInput[i] = 0; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { - m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE); - m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE); - m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ]; - m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ]; - for (int i = 0; i < m_nFilterPoints; i++) - m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0; - for (int i = 0; i < m_nOutputPoints; i++) - m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0; + m_adComplexFftInput = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints)); + m_adComplexFftOutput = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints)); + m_complexPlanForward = fftw_plan_dft_1d (m_nFilterPoints, m_adComplexFftInput, m_adComplexFftOutput, FFTW_FORWARD, FFTW_ESTIMATE); + m_adComplexFftSignal = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints)); + m_adComplexFftBackwardOutput = static_cast(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints)); + m_complexPlanBackward = fftw_plan_dft_1d (m_nOutputPoints, m_adComplexFftSignal, m_adComplexFftBackwardOutput, FFTW_BACKWARD, FFTW_ESTIMATE); + + for (i = 0; i < m_nFilterPoints; i++) + m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0; + for (i = 0; i < m_nOutputPoints; i++) + m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0; } #endif - + } ProcessSignal::~ProcessSignal (void) { - delete [] m_adFourierSinTable; - delete [] m_adFourierCosTable; - delete [] m_adFilter; + delete [] m_adFourierSinTable; + delete [] m_adFourierCosTable; + delete [] m_adFilter; #if HAVE_FFTW - if (m_idFilterMethod == FILTER_METHOD_FFTW) { - fftw_destroy_plan(m_complexPlanForward); - fftw_destroy_plan(m_complexPlanBackward); - delete [] m_adComplexFftInput; - delete [] m_adComplexFftSignal; - } - if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - rfftw_destroy_plan(m_realPlanForward); - rfftw_destroy_plan(m_realPlanBackward); - delete [] m_adRealFftInput; - delete [] m_adRealFftSignal; - } + if (m_idFilterMethod == FILTER_METHOD_FFTW) { + fftw_destroy_plan(m_complexPlanForward); + fftw_destroy_plan(m_complexPlanBackward); + fftw_free (m_adComplexFftInput); + fftw_free (m_adComplexFftOutput); + fftw_free (m_adComplexFftSignal); + fftw_free (m_adComplexFftBackwardOutput); + } + if (m_idFilterMethod == FILTER_METHOD_RFFTW) { + fftw_destroy_plan(m_realPlanForward); + fftw_destroy_plan(m_realPlanBackward); + fftw_free (m_adRealFftInput); + fftw_free (m_adRealFftOutput); + fftw_free (m_adRealFftSignal); + fftw_free (m_adRealFftBackwardOutput); + } #endif } @@ -400,12 +473,12 @@ ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName) int fmID = FILTER_METHOD_INVALID; for (int i = 0; i < s_iFilterMethodCount; i++) - if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { + if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { fmID = i; break; } - return (fmID); + return (fmID); } const char * @@ -414,7 +487,7 @@ ProcessSignal::convertFilterMethodIDToName (const int fmID) static const char *name = ""; if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterMethodName [fmID]); + return (s_aszFilterMethodName [fmID]); return (name); } @@ -425,7 +498,7 @@ ProcessSignal::convertFilterMethodIDToTitle (const int fmID) static const char *title = ""; if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterMethodTitle [fmID]); + return (s_aszFilterMethodTitle [fmID]); return (title); } @@ -437,12 +510,12 @@ ProcessSignal::convertFilterGenerationNameToID (const char* const fgName) int fgID = FILTER_GENERATION_INVALID; for (int i = 0; i < s_iFilterGenerationCount; i++) - if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) { + if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) { fgID = i; break; } - return (fgID); + return (fgID); } const char * @@ -451,7 +524,7 @@ ProcessSignal::convertFilterGenerationIDToName (const int fgID) static const char *name = ""; if (fgID >= 0 && fgID < s_iFilterGenerationCount) - return (s_aszFilterGenerationName [fgID]); + return (s_aszFilterGenerationName [fgID]); return (name); } @@ -462,108 +535,127 @@ ProcessSignal::convertFilterGenerationIDToTitle (const int fgID) static const char *name = ""; if (fgID >= 0 && fgID < s_iFilterGenerationCount) - return (s_aszFilterGenerationTitle [fgID]); + return (s_aszFilterGenerationTitle [fgID]); return (name); } void -ProcessSignal::filterSignal (const float input[], double output[]) const +ProcessSignal::filterSignal (const float constInput[], double output[]) const { + double* input = new double [m_nSignalPoints]; + int i; + for (i = 0; i < m_nSignalPoints; i++) + input[i] = constInput[i]; + + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { + for (int i = 0; i < m_nSignalPoints; i++) { + int iDetFromCenter = i - (m_nSignalPoints / 2); + input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc); + } + } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { + for (int i = 0; i < m_nSignalPoints; i++) { + int iDetFromCenter = i - (m_nSignalPoints / 2); + input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc); + } + } if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { - for (int i = 0; i < m_nSignalPoints; i++) + for (i = 0; i < m_nSignalPoints; i++) output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints); } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) { - double inputSignal[m_nFilterPoints]; - for (int i = 0; i < m_nSignalPoints; i++) + double* inputSignal = new double [m_nFilterPoints]; + for (i = 0; i < m_nSignalPoints; i++) inputSignal[i] = input[i]; - for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) + for (i = m_nSignalPoints; i < m_nFilterPoints; i++) inputSignal[i] = 0; // zeropad - complex fftSignal[m_nFilterPoints]; - finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); - for (int i = 0; i < m_nFilterPoints; i++) + std::complex* fftSignal = new std::complex [m_nFilterPoints]; + finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD); + delete inputSignal; + for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; - double inverseFourier[m_nFilterPoints]; - finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1); - for (int i = 0; i < m_nSignalPoints; i++) + double* inverseFourier = new double [m_nFilterPoints]; + finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, BACKWARD); + delete fftSignal; + for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; + delete inverseFourier; } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { - double inputSignal[m_nFilterPoints]; - for (int i = 0; i < m_nSignalPoints; i++) + double* inputSignal = new double [m_nFilterPoints]; + for (i = 0; i < m_nSignalPoints; i++) inputSignal[i] = input[i]; - for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) + for (i = m_nSignalPoints; i < m_nFilterPoints; i++) inputSignal[i] = 0; // zeropad - complex fftSignal[m_nFilterPoints]; - finiteFourierTransform (inputSignal, fftSignal, -1); - for (int i = 0; i < m_nFilterPoints; i++) + std::complex* fftSignal = new std::complex [m_nFilterPoints]; + finiteFourierTransform (inputSignal, fftSignal, FORWARD); + delete inputSignal; + for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; - double inverseFourier[m_nFilterPoints]; - finiteFourierTransform (fftSignal, inverseFourier, 1); - for (int i = 0; i < m_nSignalPoints; i++) + double* inverseFourier = new double [m_nFilterPoints]; + finiteFourierTransform (fftSignal, inverseFourier, BACKWARD); + delete fftSignal; + for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; + delete inverseFourier; } #if HAVE_FFTW else if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - for (int i = 0; i < m_nSignalPoints; i++) - m_adRealFftInput[i] = input[i]; - - fftw_real fftOutput [ m_nFilterPoints ]; - rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput); - for (int i = 0; i < m_nFilterPoints; i++) - m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i]; - for (int i = m_nFilterPoints; i < m_nOutputPoints; i++) - m_adRealFftSignal[i] = 0; - - fftw_real ifftOutput [ m_nOutputPoints ]; - rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput); - for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i]; + for (i = 0; i < m_nSignalPoints; i++) + m_adRealFftInput[i] = input[i]; + + fftw_execute (m_realPlanForward); + for (i = 0; i < m_nFilterPoints; i++) + m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i]; + for (i = m_nFilterPoints; i < m_nOutputPoints; i++) + m_adRealFftSignal[i] = 0; + + fftw_execute (m_realPlanBackward); + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = m_adRealFftBackwardOutput[i]; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { - for (int i = 0; i < m_nSignalPoints; i++) - m_adComplexFftInput[i].re = input[i]; - - fftw_complex fftOutput [ m_nFilterPoints ]; - fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput); - for (int i = 0; i < m_nFilterPoints; i++) { - m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re; - m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im; - } - fftw_complex ifftOutput [ m_nOutputPoints ]; - fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput); - for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i].re; + for (i = 0; i < m_nSignalPoints; i++) + m_adComplexFftInput[i][0] = input[i]; + + fftw_execute (m_complexPlanForward); + for (i = 0; i < m_nFilterPoints; i++) { + m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0]; + m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1]; + } + fftw_execute (m_complexPlanBackward); + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = m_adComplexFftBackwardOutput[i][0]; } #endif + delete input; } /* NAME - * convolve Discrete convolution of two functions - * - * SYNOPSIS - * r = convolve (f1, f2, dx, n, np, func_type) - * double r Convolved result - * double f1[], f2[] Functions to be convolved - * double dx Difference between successive x values - * int n Array index to center convolution about - * int np Number of points in f1 array - * int func_type EVEN or ODD or EVEN_AND_ODD function f2 - * - * NOTES - * f1 is the projection data, its indices range from 0 to np - 1. - * The index for f2, the filter, ranges from -(np-1) to (np-1). - * There are 3 ways to handle the negative vertices of f2: - * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. - * All filters used in reconstruction are even. - * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] - * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 - * for negative indices. Since f2 must range from -(np-1) to (np-1), - * if we add (np - 1) to f2's array index, then f2's index will - * range from 0 to 2 * (np - 1), and the origin, x = 0, will be - * stored at f2[np-1]. - */ - -double +* convolve Discrete convolution of two functions +* +* SYNOPSIS +* r = convolve (f1, f2, dx, n, np, func_type) +* double r Convolved result +* double f1[], f2[] Functions to be convolved +* double dx Difference between successive x values +* int n Array index to center convolution about +* int np Number of points in f1 array +* int func_type EVEN or ODD or EVEN_AND_ODD function f2 +* +* NOTES +* f1 is the projection data, its indices range from 0 to np - 1. +* The index for f2, the filter, ranges from -(np-1) to (np-1). +* There are 3 ways to handle the negative vertices of f2: +* 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. +* All filters used in reconstruction are even. +* 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] +* 3. If f2 is both ODD AND EVEN, then we must store the value of f2 +* for negative indices. Since f2 must range from -(np-1) to (np-1), +* if we add (np - 1) to f2's array index, then f2's index will +* range from 0 to 2 * (np - 1), and the origin, x = 0, will be +* stored at f2[np-1]. +*/ + +double ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const { double sum = 0.0; @@ -581,18 +673,18 @@ ProcessSignal::convolve (const double func[], const double dx, const int n, cons } -double +double ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const { double sum = 0.0; #if UNOPTIMIZED_CONVOLUTION -for (int i = 0; i < np; i++) - sum += func[i] * m_adFilter[n - i + (np - 1)]; + for (int i = 0; i < np; i++) + sum += func[i] * m_adFilter[n - i + (np - 1)]; #else -double* f2 = m_adFilter + n + (np - 1); -for (int i = 0; i < np; i++) - sum += *func++ * *f2--; + double* f2 = m_adFilter + n + (np - 1); + for (int i = 0; i < np; i++) + sum += *func++ * *f2--; #endif return (sum * dx); @@ -602,21 +694,22 @@ for (int i = 0; i < np; i++) void ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction) { - complex complexOutput[n]; + std::complex* complexOutput = new std::complex [n]; - finiteFourierTransform (input, complexOutput, n, direction); - for (int i = 0; i < n; i++) - output[i] = complexOutput[i].real(); + finiteFourierTransform (input, complexOutput, n, direction); + for (int i = 0; i < n; i++) + output[i] = complexOutput[i].real(); + delete [] complexOutput; } void -ProcessSignal::finiteFourierTransform (const double input[], complex output[], const int n, int direction) +ProcessSignal::finiteFourierTransform (const double input[], std::complex output[], const int n, int direction) { if (direction < 0) direction = -1; - else + else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { double sumReal = 0; @@ -630,25 +723,25 @@ ProcessSignal::finiteFourierTransform (const double input[], complex out sumReal /= n; sumImag /= n; } - output[i] = complex (sumReal, sumImag); + output[i] = std::complex (sumReal, sumImag); } } void -ProcessSignal::finiteFourierTransform (const complex input[], complex output[], const int n, int direction) +ProcessSignal::finiteFourierTransform (const std::complex input[], std::complex output[], const int n, int direction) { if (direction < 0) direction = -1; - else + else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { - complex sum (0,0); + std::complex sum (0,0); for (int j = 0; j < n; j++) { double angle = i * j * angleIncrement; - complex exponentTerm (cos(angle), sin(angle)); + std::complex exponentTerm (cos(angle), sin(angle)); sum += input[j] * exponentTerm; } if (direction < 0) { @@ -659,16 +752,16 @@ ProcessSignal::finiteFourierTransform (const complex input[], complex input[], double output[], const int n, int direction) +ProcessSignal::finiteFourierTransform (const std::complex input[], double output[], const int n, int direction) { if (direction < 0) direction = -1; - else + else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { - double sumReal = 0; + double sumReal = 0; for (int j = 0; j < n; j++) { double angle = i * j * angleIncrement; sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle); @@ -683,84 +776,84 @@ ProcessSignal::finiteFourierTransform (const complex input[], double out // Table-based routines void -ProcessSignal::finiteFourierTransform (const double input[], complex output[], int direction) const +ProcessSignal::finiteFourierTransform (const double input[], std::complex output[], int direction) const { if (direction < 0) direction = -1; - else + else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { double sumReal = 0, sumImag = 0; for (int j = 0; j < m_nFilterPoints; j++) { int tableIndex = i * j; if (direction > 0) { - sumReal += input[j] * m_adFourierCosTable[tableIndex]; - sumImag += input[j] * m_adFourierSinTable[tableIndex]; + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag += input[j] * m_adFourierSinTable[tableIndex]; } else { - sumReal += input[j] * m_adFourierCosTable[tableIndex]; - sumImag -= input[j] * m_adFourierSinTable[tableIndex]; + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag -= input[j] * m_adFourierSinTable[tableIndex]; } } if (direction < 0) { sumReal /= m_nFilterPoints; sumImag /= m_nFilterPoints; } - output[i] = complex (sumReal, sumImag); + output[i] = std::complex (sumReal, sumImag); } } // (a+bi) * (c + di) = (ac - bd) + (ad + bc)i void -ProcessSignal::finiteFourierTransform (const complex input[], complex output[], int direction) const +ProcessSignal::finiteFourierTransform (const std::complex input[], std::complex output[], int direction) const { if (direction < 0) direction = -1; - else + else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { double sumReal = 0, sumImag = 0; for (int j = 0; j < m_nFilterPoints; j++) { int tableIndex = i * j; if (direction > 0) { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * m_adFourierSinTable[tableIndex]; - sumImag += input[j].real() * m_adFourierSinTable[tableIndex] - + input[j].imag() * m_adFourierCosTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; } else { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * -m_adFourierSinTable[tableIndex]; - sumImag += input[j].real() * -m_adFourierSinTable[tableIndex] - + input[j].imag() * m_adFourierCosTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * -m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; } } if (direction < 0) { sumReal /= m_nFilterPoints; sumImag /= m_nFilterPoints; } - output[i] = complex (sumReal, sumImag); + output[i] = std::complex (sumReal, sumImag); } } void -ProcessSignal::finiteFourierTransform (const complex input[], double output[], int direction) const +ProcessSignal::finiteFourierTransform (const std::complex input[], double output[], int direction) const { if (direction < 0) direction = -1; - else + else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { - double sumReal = 0; + double sumReal = 0; for (int j = 0; j < m_nFilterPoints; j++) { int tableIndex = i * j; if (direction > 0) { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * m_adFourierSinTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; } else { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * -m_adFourierSinTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; } } if (direction < 0) { @@ -770,63 +863,19 @@ ProcessSignal::finiteFourierTransform (const complex input[], double out } } -// Odd Number of Points -// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2 -// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1 -// Even Number of Points -// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1) -// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1 - -void -ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) -{ - double* pdTemp = new double [n]; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[0] = pdVector[iHalfN]; - for (int i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; - for (int i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } else { // Even - int iHalfN = n / 2; - pdTemp[0] = pdVector[iHalfN]; - for (int i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + iHalfN]; - for (int i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } - - for (int i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - -void -ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n) +int +ProcessSignal::addZeropadFactor (int n, int iZeropad) { - double* pdTemp = new double [n]; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[iHalfN] = pdVector[0]; - for (int i = 0; i < iHalfN; i++) - pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; - for (int i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN + 1]; - } else { // Even - int iHalfN = n / 2; - pdTemp[iHalfN] = pdVector[0]; - for (int i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN]; - for (int i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i+1]; + if (iZeropad > 0) { + double dLogBase2 = log(n) / log(2); + int iLogBase2 = static_cast(floor (dLogBase2)); + int iPaddedN = 1 << (iLogBase2 + iZeropad); +#ifdef DEBUG + sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN); +#endif + return iPaddedN; } - for (int i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; + return n; } -