diff options
Diffstat (limited to 'gsl-1.9/fft/c_pass_n.c')
| -rw-r--r-- | gsl-1.9/fft/c_pass_n.c | 204 |
1 files changed, 204 insertions, 0 deletions
diff --git a/gsl-1.9/fft/c_pass_n.c b/gsl-1.9/fft/c_pass_n.c new file mode 100644 index 0000000..a4ec823 --- /dev/null +++ b/gsl-1.9/fft/c_pass_n.c @@ -0,0 +1,204 @@ +/* fft/c_pass_n.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. + */ + +static int +FUNCTION(fft_complex,pass_n) (BASE in[], + const size_t istride, + BASE out[], + const size_t ostride, + const gsl_fft_direction sign, + const size_t factor, + const size_t product, + const size_t n, + const TYPE(gsl_complex) twiddle[]) +{ + size_t i = 0, j = 0; + size_t k, k1; + + const size_t m = n / factor; + const size_t q = n / product; + const size_t p_1 = product / factor; + const size_t jump = (factor - 1) * p_1; + + size_t e, e1; + + for (i = 0; i < m; i++) + { + REAL(out,ostride,i) = REAL(in,istride,i); + IMAG(out,ostride,i) = IMAG(in,istride,i); + } + + for (e = 1; e < (factor - 1) / 2 + 1; e++) + { + for (i = 0; i < m; i++) + { + const size_t idx = i + e * m; + const size_t idxc = i + (factor - e) * m; + REAL(out,ostride,idx) = REAL(in,istride,idx) + REAL(in,istride,idxc); + IMAG(out,ostride,idx) = IMAG(in,istride,idx) + IMAG(in,istride,idxc); + REAL(out,ostride,idxc) = REAL(in,istride,idx) - REAL(in,istride,idxc); + IMAG(out,ostride,idxc) = IMAG(in,istride,idx) - IMAG(in,istride,idxc); + } + } + + /* e = 0 */ + + for (i=0 ; i<m; i++) + { + REAL(in,istride,i) = REAL(out,ostride,i); + IMAG(in,istride,i) = IMAG(out,ostride,i); + } + + for (e1 = 1; e1 < (factor - 1) / 2 + 1; e1++) + { + for (i = 0; i < m; i++) + { + REAL(in,istride,i) += REAL(out,ostride,i + e1*m) ; + IMAG(in,istride,i) += IMAG(out,ostride,i + e1*m) ; + } + } + + for (e = 1; e < (factor-1)/2 + 1; e++) + { + size_t idx = e*q ; + const size_t idx_step = e * q ; + ATOMIC w_real, w_imag ; + + const size_t em = e * m ; + const size_t ecm = (factor - e) * m ; + + for (i = 0; i < m; i++) + { + REAL(in,istride,i+em) = REAL(out,ostride,i) ; + IMAG(in,istride,i+em) = IMAG(out,ostride,i) ; + REAL(in,istride,i+ecm) = REAL(out,ostride,i) ; + IMAG(in,istride,i+ecm) = IMAG(out,ostride,i) ; + } + + for (e1 = 1; e1 < (factor - 1) / 2 + 1; e1++) + { + if (idx == 0) { + w_real = 1 ; + w_imag = 0 ; + } else { + if (sign == gsl_fft_forward) { + w_real = GSL_REAL(twiddle[idx - 1]) ; + w_imag = GSL_IMAG(twiddle[idx - 1]) ; + } else { + w_real = GSL_REAL(twiddle[idx - 1]) ; + w_imag = -GSL_IMAG(twiddle[idx - 1]) ; + } + } + + for (i = 0; i < m; i++) + { + const ATOMIC xp_real = REAL(out,ostride,i + e1 * m); + const ATOMIC xp_imag = IMAG(out,ostride,i + e1 * m); + const ATOMIC xm_real = REAL(out,ostride,i + (factor - e1) *m); + const ATOMIC xm_imag = IMAG(out,ostride,i + (factor - e1) *m); + + const ATOMIC ap = w_real * xp_real ; + const ATOMIC am = w_imag * xm_imag ; + + ATOMIC sum_real = ap - am; + ATOMIC sumc_real = ap + am; + + const ATOMIC bp = w_real * xp_imag ; + const ATOMIC bm = w_imag * xm_real ; + + ATOMIC sum_imag = bp + bm; + ATOMIC sumc_imag = bp - bm; + + REAL(in,istride,i + em) += sum_real; + IMAG(in,istride,i + em) += sum_imag; + REAL(in,istride,i + ecm) += sumc_real; + IMAG(in,istride,i + ecm) += sumc_imag; + } + idx += idx_step ; + idx %= factor * q ; + } + } + + i = 0; + j = 0; + + /* k = 0 */ + for (k1 = 0; k1 < p_1; k1++) + { + REAL(out,ostride,k1) = REAL(in,istride,k1); + IMAG(out,ostride,k1) = IMAG(in,istride,k1); + } + + for (e1 = 1; e1 < factor; e1++) + { + for (k1 = 0; k1 < p_1; k1++) + { + REAL(out,ostride,k1 + e1 * p_1) = REAL(in,istride,k1 + e1 * m) ; + IMAG(out,ostride,k1 + e1 * p_1) = IMAG(in,istride,k1 + e1 * m) ; + } + } + + i = p_1 ; + j = product ; + + for (k = 1; k < q; k++) + { + for (k1 = 0; k1 < p_1; k1++) + { + REAL(out,ostride,j) = REAL(in,istride,i); + IMAG(out,ostride,j) = IMAG(in,istride,i); + i++; + j++; + } + j += jump; + } + + i = p_1 ; + j = product ; + + for (k = 1; k < q; k++) + { + for (k1 = 0; k1 < p_1; k1++) + { + for (e1 = 1; e1 < factor; e1++) + { + ATOMIC x_real = REAL(in, istride,i + e1 * m); + ATOMIC x_imag = IMAG(in, istride,i + e1 * m); + + ATOMIC w_real, w_imag ; + if (sign == gsl_fft_forward) { + w_real = GSL_REAL(twiddle[(e1-1)*q + k-1]) ; + w_imag = GSL_IMAG(twiddle[(e1-1)*q + k-1]) ; + } else { + w_real = GSL_REAL(twiddle[(e1-1)*q + k-1]) ; + w_imag = -GSL_IMAG(twiddle[(e1-1)*q + k-1]) ; + } + + REAL(out,ostride,j + e1 * p_1) = w_real * x_real - w_imag * x_imag; + IMAG(out,ostride,j + e1 * p_1) = w_real * x_imag + w_imag * x_real; + } + i++; + j++; + } + j += jump; + } + + return 0; +} + |