/*****************************************************************************
** 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
#include "ct.h"
#ifdef HAVE_WXWINDOWS
-#include "dlgezplot.h"
+#include "nographics.h"
#endif
// FilterMethod ID/Names
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*);
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*);
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<int>(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);
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
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;
}
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<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
+ m_adRealFftOutput = static_cast<double*>(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<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
+ m_adRealFftBackwardOutput = static_cast<double*>(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_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
+ m_adComplexFftOutput = static_cast<fftw_complex*>(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_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
+ m_adComplexFftBackwardOutput = static_cast<fftw_complex*>(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
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
}
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;
}
}
+
+int
+ProcessSignal::addZeropadFactor (int n, int iZeropad)
+{
+ if (iZeropad > 0) {
+ double dLogBase2 = log(n) / log(2);
+ int iLogBase2 = static_cast<int>(floor (dLogBase2));
+ int iPaddedN = 1 << (iLogBase2 + iZeropad);
+#ifdef DEBUG
+ sys_error (ERR_TRACE, "Zeropadding %d to %d", n, iPaddedN);
+#endif
+ return iPaddedN;
+ }
+
+ return n;
+}
projects / ctsim.git / blobdiff