+/* 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 FunctionSymmetry func_type) const
+{
+ double sum = 0.0;
+
+ if (func_type == FUNC_BOTH) {
+#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
+ } else if (func_type == FUNC_EVEN) {
+ for (int i = 0; i < np; i++) {
+ int k = abs (n - i);
+ sum += func[i] * m_vecFilter[k];
+ }
+ } else if (func_type == FUNC_ODD) {
+ for (int i = 0; i < np; i++) {
+ int k = n - i;
+ if (k < 0)
+ sum -= func[i] * m_vecFilter[k];
+ else
+ sum += func[i] * m_vecFilter[k];
+ }
+ } else
+ sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
+
+ return (sum * dx);
+}
+
+
+double
+SignalFilter::convolve (const float func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
+{
+ double sum = 0.0;
+
+ if (func_type == FUNC_BOTH) {
+#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
+ } else if (func_type == FUNC_EVEN) {
+ for (int i = 0; i < np; i++) {
+ int k = abs (n - i);
+ sum += func[i] * m_vecFilter[k];
+ }
+ } else if (func_type == FUNC_ODD) {
+ for (int i = 0; i < np; i++) {
+ int k = n - i;
+ if (k < 0)
+ sum -= func[i] * m_vecFilter[k];
+ else
+ sum += func[i] * m_vecFilter[k];
+ }
+ } else
+ sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
+
+ return (sum * dx);
+}
+