summaryrefslogtreecommitdiff
path: root/gsl-1.9/fft/c_radix2.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/fft/c_radix2.c')
-rw-r--r--gsl-1.9/fft/c_radix2.c316
1 files changed, 316 insertions, 0 deletions
diff --git a/gsl-1.9/fft/c_radix2.c b/gsl-1.9/fft/c_radix2.c
new file mode 100644
index 0000000..f8dcf43
--- /dev/null
+++ b/gsl-1.9/fft/c_radix2.c
@@ -0,0 +1,316 @@
+/* fft/c_radix2.c
+ *
+ * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+int
+FUNCTION(gsl_fft_complex,radix2_forward) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride, const size_t n)
+{
+ gsl_fft_direction sign = gsl_fft_forward;
+ int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
+ return status;
+}
+
+int
+FUNCTION(gsl_fft_complex,radix2_backward) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride, const size_t n)
+{
+ gsl_fft_direction sign = gsl_fft_backward;
+ int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
+ return status;
+}
+
+int
+FUNCTION(gsl_fft_complex,radix2_inverse) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride, const size_t n)
+{
+ gsl_fft_direction sign = gsl_fft_backward;
+ int status = FUNCTION(gsl_fft_complex,radix2_transform) (data, stride, n, sign);
+
+ if (status)
+ {
+ return status;
+ }
+
+ /* normalize inverse fft with 1/n */
+
+ {
+ const ATOMIC norm = 1.0 / n;
+ size_t i;
+ for (i = 0; i < n; i++)
+ {
+ REAL(data,stride,i) *= norm;
+ IMAG(data,stride,i) *= norm;
+ }
+ }
+
+ return status;
+}
+
+
+
+int
+FUNCTION(gsl_fft_complex,radix2_transform) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride,
+ const size_t n,
+ const gsl_fft_direction sign)
+{
+ int result ;
+ size_t dual;
+ size_t bit;
+ size_t logn = 0;
+ int status;
+
+ if (n == 1) /* identity operation */
+ {
+ return 0 ;
+ }
+
+ /* make sure that n is a power of 2 */
+
+ result = fft_binary_logn(n) ;
+
+ if (result == -1)
+ {
+ GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
+ }
+ else
+ {
+ logn = result ;
+ }
+
+ /* bit reverse the ordering of input data for decimation in time algorithm */
+
+ status = FUNCTION(fft_complex,bitreverse_order) (data, stride, n, logn) ;
+
+ /* apply fft recursion */
+
+ dual = 1;
+
+ for (bit = 0; bit < logn; bit++)
+ {
+ ATOMIC w_real = 1.0;
+ ATOMIC w_imag = 0.0;
+
+ const double theta = 2.0 * ((int) sign) * M_PI / (2.0 * (double) dual);
+
+ const ATOMIC s = sin (theta);
+ const ATOMIC t = sin (theta / 2.0);
+ const ATOMIC s2 = 2.0 * t * t;
+
+ size_t a, b;
+
+ /* a = 0 */
+
+ for (b = 0; b < n; b += 2 * dual)
+ {
+ const size_t i = b ;
+ const size_t j = b + dual;
+
+ const ATOMIC z1_real = REAL(data,stride,j) ;
+ const ATOMIC z1_imag = IMAG(data,stride,j) ;
+
+ const ATOMIC wd_real = z1_real ;
+ const ATOMIC wd_imag = z1_imag ;
+
+ REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
+ IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
+ REAL(data,stride,i) += wd_real;
+ IMAG(data,stride,i) += wd_imag;
+ }
+
+ /* a = 1 .. (dual-1) */
+
+ for (a = 1; a < dual; a++)
+ {
+
+ /* trignometric recurrence for w-> exp(i theta) w */
+
+ {
+ const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
+ const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
+ w_real = tmp_real;
+ w_imag = tmp_imag;
+ }
+
+ for (b = 0; b < n; b += 2 * dual)
+ {
+ const size_t i = b + a;
+ const size_t j = b + a + dual;
+
+ const ATOMIC z1_real = REAL(data,stride,j) ;
+ const ATOMIC z1_imag = IMAG(data,stride,j) ;
+
+ const ATOMIC wd_real = w_real * z1_real - w_imag * z1_imag;
+ const ATOMIC wd_imag = w_real * z1_imag + w_imag * z1_real;
+
+ REAL(data,stride,j) = REAL(data,stride,i) - wd_real;
+ IMAG(data,stride,j) = IMAG(data,stride,i) - wd_imag;
+ REAL(data,stride,i) += wd_real;
+ IMAG(data,stride,i) += wd_imag;
+ }
+ }
+ dual *= 2;
+ }
+
+ return 0;
+
+}
+
+
+int
+FUNCTION(gsl_fft_complex,radix2_dif_forward) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride,
+ const size_t n)
+{
+ gsl_fft_direction sign = gsl_fft_forward;
+ int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
+ return status;
+}
+
+int
+FUNCTION(gsl_fft_complex,radix2_dif_backward) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride,
+ const size_t n)
+{
+ gsl_fft_direction sign = gsl_fft_backward;
+ int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
+ return status;
+}
+
+int
+FUNCTION(gsl_fft_complex,radix2_dif_inverse) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride,
+ const size_t n)
+{
+ gsl_fft_direction sign = gsl_fft_backward;
+ int status = FUNCTION(gsl_fft_complex,radix2_dif_transform) (data, stride, n, sign);
+
+ if (status)
+ {
+ return status;
+ }
+
+ /* normalize inverse fft with 1/n */
+
+ {
+ const ATOMIC norm = 1.0 / n;
+ size_t i;
+ for (i = 0; i < n; i++)
+ {
+ REAL(data,stride,i) *= norm;
+ IMAG(data,stride,i) *= norm;
+ }
+ }
+
+ return status;
+}
+
+int
+FUNCTION(gsl_fft_complex,radix2_dif_transform) (TYPE(gsl_complex_packed_array) data,
+ const size_t stride,
+ const size_t n,
+ const gsl_fft_direction sign)
+{
+ int result ;
+ size_t dual;
+ size_t bit;
+ size_t logn = 0;
+ int status;
+
+ if (n == 1) /* identity operation */
+ {
+ return 0 ;
+ }
+
+ /* make sure that n is a power of 2 */
+
+ result = fft_binary_logn(n) ;
+
+ if (result == -1)
+ {
+ GSL_ERROR ("n is not a power of 2", GSL_EINVAL);
+ }
+ else
+ {
+ logn = result ;
+ }
+
+ /* apply fft recursion */
+
+ dual = n / 2;
+
+ for (bit = 0; bit < logn; bit++)
+ {
+ ATOMIC w_real = 1.0;
+ ATOMIC w_imag = 0.0;
+
+ const double theta = 2.0 * ((int) sign) * M_PI / ((double) (2 * dual));
+
+ const ATOMIC s = sin (theta);
+ const ATOMIC t = sin (theta / 2.0);
+ const ATOMIC s2 = 2.0 * t * t;
+
+ size_t a, b;
+
+ for (b = 0; b < dual; b++)
+ {
+ for (a = 0; a < n; a+= 2 * dual)
+ {
+ const size_t i = b + a;
+ const size_t j = b + a + dual;
+
+ const ATOMIC t1_real = REAL(data,stride,i) + REAL(data,stride,j);
+ const ATOMIC t1_imag = IMAG(data,stride,i) + IMAG(data,stride,j);
+ const ATOMIC t2_real = REAL(data,stride,i) - REAL(data,stride,j);
+ const ATOMIC t2_imag = IMAG(data,stride,i) - IMAG(data,stride,j);
+
+ REAL(data,stride,i) = t1_real;
+ IMAG(data,stride,i) = t1_imag;
+ REAL(data,stride,j) = w_real*t2_real - w_imag * t2_imag;
+ IMAG(data,stride,j) = w_real*t2_imag + w_imag * t2_real;
+ }
+
+ /* trignometric recurrence for w-> exp(i theta) w */
+
+ {
+ const ATOMIC tmp_real = w_real - s * w_imag - s2 * w_real;
+ const ATOMIC tmp_imag = w_imag + s * w_real - s2 * w_imag;
+ w_real = tmp_real;
+ w_imag = tmp_imag;
+ }
+ }
+ dual /= 2;
+ }
+
+ /* bit reverse the ordering of output data for decimation in
+ frequency algorithm */
+
+ status = FUNCTION(fft_complex,bitreverse_order)(data, stride, n, logn) ;
+
+ return 0;
+
+}
+
+
+
+
+
+
+
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 10:33:13 +0000