X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprocsignal.cpp;h=2206ae1838cf41ebd3f2e13eb34353fdb060670e;hp=c503096945d4dc48ca997eb32a1cd9fc818cb530;hb=d42d3d062dd1aca92b5a2552a1f474aab0bee610;hpb=c953cbb6ffc2fd50e736230f4e6976a025983cff diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index c503096..2206ae1 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -1,15 +1,15 @@ /***************************************************************************** ** File IDENTIFICATION ** -** Name: filter.cpp -** Purpose: Routines for signal-procesing filters -** Progammer: Kevin Rosenberg -** Date Started: Aug 1984 +** Name: procsignal.cpp +** Purpose: Routines for processing signals and projections +** Progammer: Kevin Rosenberg +** Date Started: Aug 1984 ** ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: procsignal.cpp,v 1.27 2001/03/01 07:30:49 kevin Exp $ +** $Id$ ** ** 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 @@ -28,7 +28,7 @@ #include "ct.h" #ifdef HAVE_WXWINDOWS -#include "dlgezplot.h" +#include "nographics.h" #endif // FilterMethod ID/Names @@ -42,23 +42,23 @@ const int ProcessSignal::FILTER_METHOD_FFTW = 4; const int ProcessSignal::FILTER_METHOD_RFFTW =5 ; #endif const char* const ProcessSignal::s_aszFilterMethodName[] = { - {"convolution"}, - {"fourier"}, - {"fouier-table"}, - {"fft"}, + "convolution", + "fourier", + "fouier-table", + "fft", #if HAVE_FFTW - {"fftw"}, - {"rfftw"}, + "fftw", + "rfftw", #endif }; const char* const ProcessSignal::s_aszFilterMethodTitle[] = { - {"Convolution"}, - {"Fourier"}, - {"Fouier Trigometric Table"}, - {"FFT"}, + "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*); @@ -68,12 +68,12 @@ const int ProcessSignal::FILTER_GENERATION_INVALID = -1; const int ProcessSignal::FILTER_GENERATION_DIRECT = 0; const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1; const char* const ProcessSignal::s_aszFilterGenerationName[] = { - {"direct"}, - {"inverse-fourier"}, + "direct", + "inverse-fourier", }; const char* const ProcessSignal::s_aszFilterGenerationTitle[] = { - {"Direct"}, - {"Inverse Fourier"}, + "Direct", + "Inverse Fourier", }; const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*); @@ -258,22 +258,10 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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 defined(DEBUG) || defined(_DEBUG) - if (m_traceLevel >= Trace::TRACE_CONSOLE) - sys_error (ERR_TRACE, "nFilterPoints = %d", m_nFilterPoints); -#endif - } + m_nFilterPoints = addZeropadFactor (m_nSignalPoints, m_iZeropad); m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - if (m_nFilterPoints % 2) { // Odd + 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); @@ -297,7 +285,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw dlgEZPlot.ShowModal(); } #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 @@ -329,9 +317,9 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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; @@ -430,21 +418,26 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw } 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 ]; + 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 ]; + 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].re = m_adComplexFftInput[i].im = 0; + m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0; for (i = 0; i < m_nOutputPoints; i++) - m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0; + m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0; } #endif @@ -460,14 +453,18 @@ ProcessSignal::~ProcessSignal (void) if (m_idFilterMethod == FILTER_METHOD_FFTW) { fftw_destroy_plan(m_complexPlanForward); fftw_destroy_plan(m_complexPlanBackward); - delete [] m_adComplexFftInput; - delete [] m_adComplexFftSignal; + fftw_free (m_adComplexFftInput); + fftw_free (m_adComplexFftOutput); + fftw_free (m_adComplexFftSignal); + fftw_free (m_adComplexFftBackwardOutput); } if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - rfftw_destroy_plan(m_realPlanForward); - rfftw_destroy_plan(m_realPlanBackward); - delete [] m_adRealFftInput; - delete [] m_adRealFftSignal; + 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 } @@ -607,35 +604,27 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const for (i = 0; i < m_nSignalPoints; i++) m_adRealFftInput[i] = input[i]; - fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ]; - rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput); + fftw_execute (m_realPlanForward); for (i = 0; i < m_nFilterPoints; i++) - m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i]; - delete [] fftOutput; + m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i]; for (i = m_nFilterPoints; i < m_nOutputPoints; i++) m_adRealFftSignal[i] = 0; - fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ]; - rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput); + fftw_execute (m_realPlanBackward); for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i]; - delete [] ifftOutput; + output[i] = m_adRealFftBackwardOutput[i]; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { for (i = 0; i < m_nSignalPoints; i++) - m_adComplexFftInput[i].re = input[i]; + m_adComplexFftInput[i][0] = input[i]; - fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ]; - fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput); + fftw_execute (m_complexPlanForward); for (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; + m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0]; + m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1]; } - delete [] fftOutput; - fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ]; - fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput); + fftw_execute (m_complexPlanBackward); for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i].re; - delete [] ifftOutput; + output[i] = m_adComplexFftBackwardOutput[i][0]; } #endif delete input; @@ -876,3 +865,19 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl } } + +int +ProcessSignal::addZeropadFactor (int n, int iZeropad) +{ + 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; + } + + return n; +}