/*****************************************************************************
** 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.29 2001/03/18 18:08:25 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_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;
if (iZeropad > 0) {
double dLogBase2 = log(n) / log(2);
int iLogBase2 = static_cast<int>(floor (dLogBase2));
- if (dLogBase2 != floor(dLogBase2))
- iLogBase2++; // raise up to next power of 2
- n = 1 << (iLogBase2 + (iZeropad - 1));
+ 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