-/* NAME
- * convolve Discrete convolution of two functions
- *
- * SYNOPSIS
- * r = convolve (f1, f2, dx, n, np, func_type)
- * double r Convolved result
- * double f1[], f2[] Functions to be convolved
- * double dx Difference between successive x values
- * int n Array index to center convolution about
- * int np Number of points in f1 array
- * int func_type EVEN or ODD or EVEN_AND_ODD function f2
- *
- * NOTES
- * f1 is the projection data, its indices range from 0 to np - 1.
- * The index for f2, the filter, ranges from -(np-1) to (np-1).
- * There are 3 ways to handle the negative vertices of f2:
- * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n].
- * All filters used in reconstruction are even.
- * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n]
- * 3. If f2 is both ODD AND EVEN, then we must store the value of f2
- * for negative indices. Since f2 must range from -(np-1) to (np-1),
- * if we add (np - 1) to f2's array index, then f2's index will
- * range from 0 to 2 * (np - 1), and the origin, x = 0, will be
- * stored at f2[np-1].
- */
-
-double
-SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const
-{
- double sum = 0.0;
-
-#if UNOPTIMIZED_CONVOLUTION
- for (int i = 0; i < np; i++)
- sum += func[i] * m_vecFilter[n - i + (np - 1)];
-#else
- double* f2 = m_vecFilter + n + (np - 1);
- for (int i = 0; i < np; i++)
- sum += *func++ * *f2--;
-#endif
-
- return (sum * dx);
-}
-
-
-double
-SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const
-{
- double sum = 0.0;
-
-#if UNOPTIMIZED_CONVOLUTION
-for (int i = 0; i < np; i++)
- sum += func[i] * m_vecFilter[n - i + (np - 1)];
-#else
-double* f2 = m_vecFilter + n + (np - 1);
-for (int i = 0; i < np; i++)
- sum += *func++ * *f2--;
-#endif
-
- return (sum * dx);
-}
-
-
-void
-SignalFilter::finiteFourierTransform (const float input[], complex<double> output[], const int n, int direction)
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- double angleIncrement = 2 * PI / n;
- for (int i = 0; i < n; i++) {
- double sumReal = 0;
- double sumImag = 0;
- for (int j = 0; j < n; j++) {
- double angle = i * j * angleIncrement * direction;
- sumReal += input[i] * cos(angle);
- sumImag += input[i] * sin(angle);
- }
- if (direction > 0) {
- sumReal /= n;
- sumImag /= n;
- }
- output[i] = complex<double> (sumReal, sumImag);
- }
-}
-
-void
-SignalFilter::finiteFourierTransform (const float input[], complex<double> output[], int direction) const
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- double angleIncrement = 2 * PI / m_nSignalPoints;
- for (int i = 0; i < m_nSignalPoints; i++) {
- double sumReal = 0, sumImag = 0;
- for (int j = 0; j < m_nSignalPoints; j++) {
- int tableIndex = i * j;
- if (direction > 0) {
- sumReal += input[i] * m_vecFourierCosTable[tableIndex];
- sumImag += input[i] * m_vecFourierSinTable[tableIndex];
- } else {
- sumReal += input[i] * m_vecFourierCosTable[tableIndex];
- sumImag -= input[i] * m_vecFourierSinTable[tableIndex];
- }
- }
- if (direction > 0) {
- sumReal /= m_nSignalPoints;
- sumImag /= m_nSignalPoints;
- }
- output[i] = complex<double> (sumReal, sumImag);
- }
-}