-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 double input[], complex<double> output[], const int n, int direction)
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- double angleIncrement = direction * 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;
- sumReal += input[j] * cos(angle);
- sumImag += input[j] * sin(angle);
- }
- if (direction < 0) {
- sumReal /= n;
- sumImag /= n;
- }
- output[i] = complex<double> (sumReal, sumImag);
- }
-}
-
-
-void
-SignalFilter::finiteFourierTransform (const complex<double> input[], complex<double> output[], const int n, int direction)
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- double angleIncrement = direction * 2 * PI / n;
- for (int i = 0; i < n; i++) {
- complex<double> sum (0,0);
- for (int j = 0; j < n; j++) {
- double angle = i * j * angleIncrement;
- complex<double> exponentTerm (cos(angle), sin(angle));
- sum += input[j] * exponentTerm;
- }
- if (direction < 0) {
- sum /= n;
- }
- output[i] = sum;
- }
-}
-
-void
-SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], const int n, int direction)
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- double angleIncrement = direction * 2 * PI / n;
- for (int i = 0; i < n; i++) {
- double sumReal = 0;
- for (int j = 0; j < n; j++) {
- double angle = i * j * angleIncrement;
- sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle);
- }
- if (direction < 0) {
- sumReal /= n;
- }
- output[i] = sumReal;
- }
-}
-
-void
-SignalFilter::finiteFourierTransform (const double input[], complex<double> output[], int direction) const
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- for (int i = 0; i < m_nFilterPoints; i++) {
- double sumReal = 0, sumImag = 0;
- for (int j = 0; j < m_nFilterPoints; j++) {
- int tableIndex = i * j;
- if (direction > 0) {
- sumReal += input[j] * m_vecFourierCosTable[tableIndex];
- sumImag += input[j] * m_vecFourierSinTable[tableIndex];
- } else {
- sumReal += input[j] * m_vecFourierCosTable[tableIndex];
- sumImag -= input[j] * m_vecFourierSinTable[tableIndex];
- }
- }
- if (direction < 0) {
- sumReal /= m_nFilterPoints;
- sumImag /= m_nFilterPoints;
- }
- output[i] = complex<double> (sumReal, sumImag);
- }
-}
-
-// (a+bi) * (c + di) = (ac - bd) + (ad + bc)i
-void
-SignalFilter::finiteFourierTransform (const complex<double> input[], complex<double> output[], int direction) const
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- for (int i = 0; i < m_nFilterPoints; i++) {
- double sumReal = 0, sumImag = 0;
- for (int j = 0; j < m_nFilterPoints; j++) {
- int tableIndex = i * j;
- if (direction > 0) {
- sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
- - input[j].imag() * m_vecFourierSinTable[tableIndex];
- sumImag += input[j].real() * m_vecFourierSinTable[tableIndex]
- + input[j].imag() * m_vecFourierCosTable[tableIndex];
- } else {
- sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
- - input[j].imag() * -m_vecFourierSinTable[tableIndex];
- sumImag += input[j].real() * -m_vecFourierSinTable[tableIndex]
- + input[j].imag() * m_vecFourierCosTable[tableIndex];
- }
- }
- if (direction < 0) {
- sumReal /= m_nFilterPoints;
- sumImag /= m_nFilterPoints;
- }
- output[i] = complex<double> (sumReal, sumImag);
- }
-}
-
-void
-SignalFilter::finiteFourierTransform (const complex<double> input[], double output[], int direction) const
-{
- if (direction < 0)
- direction = -1;
- else
- direction = 1;
-
- for (int i = 0; i < m_nFilterPoints; i++) {
- double sumReal = 0;
- for (int j = 0; j < m_nFilterPoints; j++) {
- int tableIndex = i * j;
- if (direction > 0) {
- sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
- - input[j].imag() * m_vecFourierSinTable[tableIndex];
- } else {
- sumReal += input[j].real() * m_vecFourierCosTable[tableIndex]
- - input[j].imag() * -m_vecFourierSinTable[tableIndex];
- }
- }
- if (direction < 0) {
- sumReal /= m_nFilterPoints;
- }
- output[i] = sumReal;
- }
-}