** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: procsignal.cpp,v 1.1 2000/08/19 23:00:05 kevin Exp $
+** $Id: procsignal.cpp,v 1.3 2000/08/25 15:59:13 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
};
const char* ProcessSignal::s_aszFilterMethodTitle[] = {
{"Convolution"},
- {"Direct Fourier"},
- {"Fouier Trigometric Table Lookout"},
+ {"Fourier"},
+ {"Fouier Trigometric Table"},
{"FFT"},
#if HAVE_FFTW
{"FFTW"},
// 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)
+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)
: m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
{
m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
return;
}
- init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor);
+ init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel);
}
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)
+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)
{
m_idFilter = idFilter;
m_idDomain = idDomain;
m_fail = true;
return;
}
- m_traceLevel = TRACE_NONE;
+ m_traceLevel = iTraceLevel;
m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration);
m_dBandwidth = dBandwidth;
if (! m_bFrequencyFiltering) {
if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
- m_nFilterPoints = 2 * m_nSignalPoints - 1;
+ 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);
m_adFilter = new double[ m_nFilterPoints ];
filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
} else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
- m_nFilterPoints = m_nSignalPoints;
+ 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];
- double adInverseFilter [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);
- ProcessSignal::finiteFourierTransform (adFrequencyFilter, adInverseFilter, m_nFilterPoints, 1);
- for (int i = 0; i < m_nFilterPoints; i++)
- m_adFilter [i] = adInverseFilter[i];
+ 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;
+ }
}
}
// Frequency-based filtering
else if (m_bFrequencyFiltering) {
- // calculate number of filter points with zeropadding
- m_nFilterPoints = m_nSignalPoints;
- if (m_iZeropad > 0) {
- double logBase2 = log(m_nSignalPoints) / log(2);
- int nextPowerOf2 = static_cast<int>(floor(logBase2));
- if (logBase2 != floor(logBase2))
- nextPowerOf2++;
- nextPowerOf2 += (m_iZeropad - 1);
- m_nFilterPoints = 1 << nextPowerOf2;
- if (m_traceLevel >= TRACE_TEXT)
+ 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 (m_traceLevel >= TRACE_TEXT)
cout << "nFilterPoints = " << m_nFilterPoints << endl;
- }
- m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
+ }
+ m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor;
- if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
+ 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);
- 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);
- shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints);
+ } 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;
+ }
+
+ 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 ();
+ }
+ 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 ();
+ }
} else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) {
- m_nFilterPoints = 2 * m_nSignalPoints - 1;
- m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1);
- m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1);
- m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1);
- double adSpatialFilter [m_nFilterPoints];
- double adInverseFilter [m_nFilterPoints];
- SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL);
- filter.copyFilterData (adSpatialFilter, 0, m_nFilterPoints);
- m_adFilter = new double [m_nFilterPoints];
- finiteFourierTransform (adSpatialFilter, adInverseFilter, m_nFilterPoints, -1);
- for (int i = 0; i < m_nFilterPoints; i++)
- m_adFilter [i] = adInverseFilter[i];
+ // calculate number of filter points with zeropadding
+ int nSpatialPoints = 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) / (nSpatialPoints - 1);
+ m_nFilterPoints = nSpatialPoints;
+ if (m_iZeropad > 0) {
+ double logBase2 = log(nSpatialPoints) / log(2);
+ int nextPowerOf2 = static_cast<int>(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);
+ 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 ();
+ }
+ for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
+ adSpatialFilter[i] = 0;
+ m_adFilter = new double [m_nFilterPoints];
+ complex<double> 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 ();
+ }
+ }
+ }
+
// 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;
{
delete [] m_adFourierSinTable;
delete [] m_adFourierCosTable;
+ delete [] m_adFilter;
#if HAVE_FFTW
if (m_idFilterMethod == FILTER_METHOD_FFTW) {
finiteFourierTransform (input, complexOutput, n, direction);
for (int i = 0; i < n; i++)
- output[i] = abs(complexOutput[n]);
+ output[i] = complexOutput[i].real();
}
void
int iHalfN = (n - 1) / 2;
pdTemp[0] = pdVector[iHalfN];
- for (int i = 1; i <= iHalfN; i++)
- pdTemp[i] = pdVector[i+iHalfN];
- for (int i = iHalfN+1; i < n; i++)
- pdTemp[i] = pdVector[i-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++)
int iHalfN = (n - 1) / 2;
pdTemp[iHalfN] = pdVector[0];
- for (int i = 1; i <= iHalfN; i++)
- pdTemp[i] = pdVector[i+iHalfN];
- for (int i = iHalfN+1; i < n; i++)
- pdTemp[i] = pdVector[i-iHalfN];
+ 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];
}
for (int i = 0; i < n; i++)
projects / ctsim.git / blobdiff