** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: procsignal.cpp,v 1.7 2000/09/07 14:29:05 kevin Exp $
+** $Id: procsignal.cpp,v 1.8 2000/12/06 01:46:43 kevin Exp $
**
** 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
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];
+ double* adFrequencyFilter = new double [m_nFilterPoints];
filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints);
#ifdef HAVE_SGP
EZPlot* pEZPlot = NULL;
}
#endif
ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1);
+ delete adFrequencyFilter;\r
#ifdef HAVE_SGP
if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order");
if (m_traceLevel >= Trace::TRACE_CONSOLE)
cout << "nFilterPoints = " << m_nFilterPoints << endl;
#endif
- double adSpatialFilter [m_nFilterPoints];
+ double* adSpatialFilter = new double [m_nFilterPoints];\r
SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints);
#ifdef HAVE_SGP
double dScale = 0.5 * sinScale * sinScale;
adSpatialFilter[i] *= dScale;
}
- }
- for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
+ }\r
+ int i;
+ for (i = nSpatialPoints; i < m_nFilterPoints; i++)
adSpatialFilter[i] = 0;
m_adFilter = new double [m_nFilterPoints];
- complex<double> acInverseFilter [m_nFilterPoints];
+ complex<double>* acInverseFilter = new complex<double> [m_nFilterPoints];\r
finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1);
- for (int i = 0; i < m_nFilterPoints; i++)
- m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
+ delete adSpatialFilter;\r
+ for (i = 0; i < m_nFilterPoints; i++)
+ m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc;
+ delete acInverseFilter;\r
#ifdef HAVE_SGP
if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) {
pEZPlot->ezset ("title Spatial Filter: Inverse");
pEZPlot->ezset ("yporigin 0.50");
pEZPlot->addCurve (m_adFilter, m_nFilterPoints);
pEZPlot->plot();
- delete pEZPlot;
+ delete pEZPlot;\r
}
#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 ];
void
ProcessSignal::filterSignal (const float constInput[], double output[]) const
{
- double input [m_nSignalPoints];
- for (int i = 0; i < m_nSignalPoints; i++)
+ double* input = new double [m_nSignalPoints];
+ int i;\r
+ for (i = 0; i < m_nSignalPoints; i++)
input[i] = constInput[i];
if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
int iDetFromCenter = i - (m_nSignalPoints / 2);
input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
}
- }
+ }\r
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<double> fftSignal[m_nFilterPoints];
- finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);
- for (int i = 0; i < m_nFilterPoints; i++)
+ complex<double>* fftSignal = new complex<double> [m_nFilterPoints];
+ finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1);\r
+ 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++)
- output[i] = inverseFourier[i];
+ double* inverseFourier = new double [m_nFilterPoints];
+ finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1);\r
+ delete fftSignal;
+ for (i = 0; i < m_nSignalPoints; i++)
+ output[i] = inverseFourier[i];\r
+ 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<double> fftSignal[m_nFilterPoints];
- finiteFourierTransform (inputSignal, fftSignal, -1);
- for (int i = 0; i < m_nFilterPoints; i++)
+ complex<double>* fftSignal = new complex<double> [m_nFilterPoints];
+ finiteFourierTransform (inputSignal, fftSignal, -1);\r
+ 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++)
- output[i] = inverseFourier[i];
+ double* inverseFourier = new double [m_nFilterPoints];
+ finiteFourierTransform (fftSignal, inverseFourier, 1);\r
+ delete fftSignal;
+ for (i = 0; i < m_nSignalPoints; i++)
+ output[i] = inverseFourier[i];\r
+ delete inverseFourier;
}
#if HAVE_FFTW
else if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
output[i] = ifftOutput[i].re;
}
-#endif
+#endif\r
+ delete input;
}
void
ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction)
{
- complex<double> complexOutput[n];
+ complex<double>* complexOutput = new complex<double> [n];
finiteFourierTransform (input, complexOutput, n, direction);
for (int i = 0; i < n; i++)
- output[i] = complexOutput[i].real();
+ output[i] = complexOutput[i].real();\r
+ delete [] complexOutput;
}
void
void
ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n)
{
- double* pdTemp = new double [n];
+ double* pdTemp = new double [n];\r
+ int i;
if (n % 2) { // Odd
int iHalfN = (n - 1) / 2;
pdTemp[0] = pdVector[iHalfN];
- for (int i = 0; i < iHalfN; i++)
+ for (i = 0; i < iHalfN; i++)
pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
- for (int i = 0; i < iHalfN; i++)
+ for (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++)
+ for (i = 0; i < iHalfN; i++)
pdTemp[i + 1] = pdVector[i + iHalfN];
- for (int i = 0; i < iHalfN - 1; i++)
+ for (i = 0; i < iHalfN - 1; i++)
pdTemp[i + iHalfN + 1] = pdVector[i];
}
- for (int i = 0; i < n; i++)
+ for (i = 0; i < n; i++)
pdVector[i] = pdTemp[i];
delete pdTemp;
}
void
ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n)
{
- double* pdTemp = new double [n];
+ double* pdTemp = new double [n];\r
+ int i;
if (n % 2) { // Odd
int iHalfN = (n - 1) / 2;
- pdTemp[iHalfN] = pdVector[0];
- for (int i = 0; i < iHalfN; i++)
+ pdTemp[iHalfN] = pdVector[0];\r
+ for (i = 0; i < iHalfN; i++)
pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
- for (int i = 0; i < iHalfN; i++)
+ for (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++)
+ for (i = 0; i < iHalfN; i++)
pdTemp[i] = pdVector[i + iHalfN];
- for (int i = 0; i < iHalfN - 1; i++)
+ for (i = 0; i < iHalfN - 1; i++)
pdTemp[i + iHalfN + 1] = pdVector[i+1];
}
- for (int i = 0; i < n; i++)
+ for (i = 0; i < n; i++)
pdVector[i] = pdTemp[i];
delete pdTemp;
}