| /* Copyright (c) 2013 Julien Pommier ( pommier@modartt.com ) |
| Copyright (c) 2020 Hayati Ayguen ( h_ayguen@web.de ) |
| |
| Based on original fortran 77 code from FFTPACKv4 from NETLIB |
| (http://www.netlib.org/fftpack), authored by Dr Paul Swarztrauber |
| of NCAR, in 1985. |
| |
| As confirmed by the NCAR fftpack software curators, the following |
| FFTPACKv5 license applies to FFTPACKv4 sources. My changes are |
| released under the same terms. |
| |
| FFTPACK license: |
| |
| http://www.cisl.ucar.edu/css/software/fftpack5/ftpk.html |
| |
| Copyright (c) 2004 the University Corporation for Atmospheric |
| Research ("UCAR"). All rights reserved. Developed by NCAR's |
| Computational and Information Systems Laboratory, UCAR, |
| www.cisl.ucar.edu. |
| |
| Redistribution and use of the Software in source and binary forms, |
| with or without modification, is permitted provided that the |
| following conditions are met: |
| |
| - Neither the names of NCAR's Computational and Information Systems |
| Laboratory, the University Corporation for Atmospheric Research, |
| nor the names of its sponsors or contributors may be used to |
| endorse or promote products derived from this Software without |
| specific prior written permission. |
| |
| - Redistributions of source code must retain the above copyright |
| notices, this list of conditions, and the disclaimer below. |
| |
| - Redistributions in binary form must reproduce the above copyright |
| notice, this list of conditions, and the disclaimer below in the |
| documentation and/or other materials provided with the |
| distribution. |
| |
| THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, |
| EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF |
| MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND |
| NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT |
| HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL, |
| EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN |
| ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN |
| CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE |
| SOFTWARE. |
| |
| |
| PFFFT : a Pretty Fast FFT. |
| |
| This file is largerly based on the original FFTPACK implementation, modified in |
| order to take advantage of SIMD instructions of modern CPUs. |
| */ |
| |
| /* |
| ChangeLog: |
| - 2011/10/02, version 1: This is the very first release of this file. |
| */ |
| |
| #include "pffft.h" |
| |
| /* detect compiler flavour */ |
| #if defined(_MSC_VER) |
| # define COMPILER_MSVC |
| #elif defined(__GNUC__) |
| # define COMPILER_GCC |
| #endif |
| |
| #include <stdlib.h> |
| #include <stdint.h> |
| #include <stdio.h> |
| #include <math.h> |
| #include <assert.h> |
| |
| #if defined(COMPILER_GCC) |
| # define ALWAYS_INLINE(return_type) inline return_type __attribute__ ((always_inline)) |
| # define NEVER_INLINE(return_type) return_type __attribute__ ((noinline)) |
| # define RESTRICT __restrict |
| # define VLA_ARRAY_ON_STACK(type__, varname__, size__) type__ varname__[size__]; |
| #elif defined(COMPILER_MSVC) |
| # define ALWAYS_INLINE(return_type) __forceinline return_type |
| # define NEVER_INLINE(return_type) __declspec(noinline) return_type |
| # define RESTRICT __restrict |
| # define VLA_ARRAY_ON_STACK(type__, varname__, size__) type__ *varname__ = (type__*)_alloca(size__ * sizeof(type__)) |
| #endif |
| |
| |
| #ifdef COMPILER_MSVC |
| #pragma warning( disable : 4244 4305 4204 4456 ) |
| #endif |
| |
| /* |
| vector support macros: the rest of the code is independant of |
| SSE/Altivec/NEON -- adding support for other platforms with 4-element |
| vectors should be limited to these macros |
| */ |
| |
| |
| // define PFFFT_SIMD_DISABLE if you want to use scalar code instead of simd code |
| //#define PFFFT_SIMD_DISABLE |
| |
| #ifdef PFFFT_FLOAT |
| #define float PFFFT_FLOAT |
| #endif |
| |
| |
| #include "pf_float.h" |
| |
| /* detect bugs with the vector support macros */ |
| void validate_pffft_simd() { |
| #ifndef PFFFT_SIMD_DISABLE |
| Vvalidate_simd(); |
| #endif |
| } |
| |
| /* SSE and co like 16-bytes aligned pointers */ |
| #define MALLOC_V4SF_ALIGNMENT 64 // with a 64-byte alignment, we are even aligned on L2 cache lines... |
| void *pffft_aligned_malloc(size_t nb_bytes) { |
| return Valigned_malloc(nb_bytes); |
| } |
| |
| void pffft_aligned_free(void *p) { |
| Valigned_free(p); |
| } |
| |
| int pffft_simd_size() { return SIMD_SZ; } |
| |
| /* |
| passf2 and passb2 has been merged here, fsign = -1 for passf2, +1 for passb2 |
| */ |
| static NEVER_INLINE(void) passf2_ps(int ido, int l1, const v4sf *cc, v4sf *ch, const float *wa1, float fsign) { |
| int k, i; |
| int l1ido = l1*ido; |
| if (ido <= 2) { |
| for (k=0; k < l1ido; k += ido, ch += ido, cc+= 2*ido) { |
| ch[0] = VADD(cc[0], cc[ido+0]); |
| ch[l1ido] = VSUB(cc[0], cc[ido+0]); |
| ch[1] = VADD(cc[1], cc[ido+1]); |
| ch[l1ido + 1] = VSUB(cc[1], cc[ido+1]); |
| } |
| } else { |
| for (k=0; k < l1ido; k += ido, ch += ido, cc += 2*ido) { |
| for (i=0; i<ido-1; i+=2) { |
| v4sf tr2 = VSUB(cc[i+0], cc[i+ido+0]); |
| v4sf ti2 = VSUB(cc[i+1], cc[i+ido+1]); |
| v4sf wr = LD_PS1(wa1[i]), wi = VMUL(LD_PS1(fsign), LD_PS1(wa1[i+1])); |
| ch[i] = VADD(cc[i+0], cc[i+ido+0]); |
| ch[i+1] = VADD(cc[i+1], cc[i+ido+1]); |
| VCPLXMUL(tr2, ti2, wr, wi); |
| ch[i+l1ido] = tr2; |
| ch[i+l1ido+1] = ti2; |
| } |
| } |
| } |
| } |
| |
| /* |
| passf3 and passb3 has been merged here, fsign = -1 for passf3, +1 for passb3 |
| */ |
| static NEVER_INLINE(void) passf3_ps(int ido, int l1, const v4sf *cc, v4sf *ch, |
| const float *wa1, const float *wa2, float fsign) { |
| static const float taur = -0.5f; |
| float taui = 0.866025403784439f*fsign; |
| int i, k; |
| v4sf tr2, ti2, cr2, ci2, cr3, ci3, dr2, di2, dr3, di3; |
| int l1ido = l1*ido; |
| float wr1, wi1, wr2, wi2; |
| assert(ido > 2); |
| for (k=0; k< l1ido; k += ido, cc+= 3*ido, ch +=ido) { |
| for (i=0; i<ido-1; i+=2) { |
| tr2 = VADD(cc[i+ido], cc[i+2*ido]); |
| cr2 = VADD(cc[i], SVMUL(taur,tr2)); |
| ch[i] = VADD(cc[i], tr2); |
| ti2 = VADD(cc[i+ido+1], cc[i+2*ido+1]); |
| ci2 = VADD(cc[i +1], SVMUL(taur,ti2)); |
| ch[i+1] = VADD(cc[i+1], ti2); |
| cr3 = SVMUL(taui, VSUB(cc[i+ido], cc[i+2*ido])); |
| ci3 = SVMUL(taui, VSUB(cc[i+ido+1], cc[i+2*ido+1])); |
| dr2 = VSUB(cr2, ci3); |
| dr3 = VADD(cr2, ci3); |
| di2 = VADD(ci2, cr3); |
| di3 = VSUB(ci2, cr3); |
| wr1=wa1[i], wi1=fsign*wa1[i+1], wr2=wa2[i], wi2=fsign*wa2[i+1]; |
| VCPLXMUL(dr2, di2, LD_PS1(wr1), LD_PS1(wi1)); |
| ch[i+l1ido] = dr2; |
| ch[i+l1ido + 1] = di2; |
| VCPLXMUL(dr3, di3, LD_PS1(wr2), LD_PS1(wi2)); |
| ch[i+2*l1ido] = dr3; |
| ch[i+2*l1ido+1] = di3; |
| } |
| } |
| } /* passf3 */ |
| |
| static NEVER_INLINE(void) passf4_ps(int ido, int l1, const v4sf *cc, v4sf *ch, |
| const float *wa1, const float *wa2, const float *wa3, float fsign) { |
| /* isign == -1 for forward transform and +1 for backward transform */ |
| |
| int i, k; |
| v4sf ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; |
| int l1ido = l1*ido; |
| if (ido == 2) { |
| for (k=0; k < l1ido; k += ido, ch += ido, cc += 4*ido) { |
| tr1 = VSUB(cc[0], cc[2*ido + 0]); |
| tr2 = VADD(cc[0], cc[2*ido + 0]); |
| ti1 = VSUB(cc[1], cc[2*ido + 1]); |
| ti2 = VADD(cc[1], cc[2*ido + 1]); |
| ti4 = VMUL(VSUB(cc[1*ido + 0], cc[3*ido + 0]), LD_PS1(fsign)); |
| tr4 = VMUL(VSUB(cc[3*ido + 1], cc[1*ido + 1]), LD_PS1(fsign)); |
| tr3 = VADD(cc[ido + 0], cc[3*ido + 0]); |
| ti3 = VADD(cc[ido + 1], cc[3*ido + 1]); |
| |
| ch[0*l1ido + 0] = VADD(tr2, tr3); |
| ch[0*l1ido + 1] = VADD(ti2, ti3); |
| ch[1*l1ido + 0] = VADD(tr1, tr4); |
| ch[1*l1ido + 1] = VADD(ti1, ti4); |
| ch[2*l1ido + 0] = VSUB(tr2, tr3); |
| ch[2*l1ido + 1] = VSUB(ti2, ti3); |
| ch[3*l1ido + 0] = VSUB(tr1, tr4); |
| ch[3*l1ido + 1] = VSUB(ti1, ti4); |
| } |
| } else { |
| for (k=0; k < l1ido; k += ido, ch+=ido, cc += 4*ido) { |
| for (i=0; i<ido-1; i+=2) { |
| float wr1, wi1, wr2, wi2, wr3, wi3; |
| tr1 = VSUB(cc[i + 0], cc[i + 2*ido + 0]); |
| tr2 = VADD(cc[i + 0], cc[i + 2*ido + 0]); |
| ti1 = VSUB(cc[i + 1], cc[i + 2*ido + 1]); |
| ti2 = VADD(cc[i + 1], cc[i + 2*ido + 1]); |
| tr4 = VMUL(VSUB(cc[i + 3*ido + 1], cc[i + 1*ido + 1]), LD_PS1(fsign)); |
| ti4 = VMUL(VSUB(cc[i + 1*ido + 0], cc[i + 3*ido + 0]), LD_PS1(fsign)); |
| tr3 = VADD(cc[i + ido + 0], cc[i + 3*ido + 0]); |
| ti3 = VADD(cc[i + ido + 1], cc[i + 3*ido + 1]); |
| |
| ch[i] = VADD(tr2, tr3); |
| cr3 = VSUB(tr2, tr3); |
| ch[i + 1] = VADD(ti2, ti3); |
| ci3 = VSUB(ti2, ti3); |
| |
| cr2 = VADD(tr1, tr4); |
| cr4 = VSUB(tr1, tr4); |
| ci2 = VADD(ti1, ti4); |
| ci4 = VSUB(ti1, ti4); |
| wr1=wa1[i], wi1=fsign*wa1[i+1]; |
| VCPLXMUL(cr2, ci2, LD_PS1(wr1), LD_PS1(wi1)); |
| wr2=wa2[i], wi2=fsign*wa2[i+1]; |
| ch[i + l1ido] = cr2; |
| ch[i + l1ido + 1] = ci2; |
| |
| VCPLXMUL(cr3, ci3, LD_PS1(wr2), LD_PS1(wi2)); |
| wr3=wa3[i], wi3=fsign*wa3[i+1]; |
| ch[i + 2*l1ido] = cr3; |
| ch[i + 2*l1ido + 1] = ci3; |
| |
| VCPLXMUL(cr4, ci4, LD_PS1(wr3), LD_PS1(wi3)); |
| ch[i + 3*l1ido] = cr4; |
| ch[i + 3*l1ido + 1] = ci4; |
| } |
| } |
| } |
| } /* passf4 */ |
| |
| /* |
| passf5 and passb5 has been merged here, fsign = -1 for passf5, +1 for passb5 |
| */ |
| static NEVER_INLINE(void) passf5_ps(int ido, int l1, const v4sf *cc, v4sf *ch, |
| const float *wa1, const float *wa2, |
| const float *wa3, const float *wa4, float fsign) { |
| static const float tr11 = .309016994374947f; |
| const float ti11 = .951056516295154f*fsign; |
| static const float tr12 = -.809016994374947f; |
| const float ti12 = .587785252292473f*fsign; |
| |
| /* Local variables */ |
| int i, k; |
| v4sf ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, |
| ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5; |
| |
| float wr1, wi1, wr2, wi2, wr3, wi3, wr4, wi4; |
| |
| #define cc_ref(a_1,a_2) cc[(a_2-1)*ido + a_1 + 1] |
| #define ch_ref(a_1,a_3) ch[(a_3-1)*l1*ido + a_1 + 1] |
| |
| assert(ido > 2); |
| for (k = 0; k < l1; ++k, cc += 5*ido, ch += ido) { |
| for (i = 0; i < ido-1; i += 2) { |
| ti5 = VSUB(cc_ref(i , 2), cc_ref(i , 5)); |
| ti2 = VADD(cc_ref(i , 2), cc_ref(i , 5)); |
| ti4 = VSUB(cc_ref(i , 3), cc_ref(i , 4)); |
| ti3 = VADD(cc_ref(i , 3), cc_ref(i , 4)); |
| tr5 = VSUB(cc_ref(i-1, 2), cc_ref(i-1, 5)); |
| tr2 = VADD(cc_ref(i-1, 2), cc_ref(i-1, 5)); |
| tr4 = VSUB(cc_ref(i-1, 3), cc_ref(i-1, 4)); |
| tr3 = VADD(cc_ref(i-1, 3), cc_ref(i-1, 4)); |
| ch_ref(i-1, 1) = VADD(cc_ref(i-1, 1), VADD(tr2, tr3)); |
| ch_ref(i , 1) = VADD(cc_ref(i , 1), VADD(ti2, ti3)); |
| cr2 = VADD(cc_ref(i-1, 1), VADD(SVMUL(tr11, tr2),SVMUL(tr12, tr3))); |
| ci2 = VADD(cc_ref(i , 1), VADD(SVMUL(tr11, ti2),SVMUL(tr12, ti3))); |
| cr3 = VADD(cc_ref(i-1, 1), VADD(SVMUL(tr12, tr2),SVMUL(tr11, tr3))); |
| ci3 = VADD(cc_ref(i , 1), VADD(SVMUL(tr12, ti2),SVMUL(tr11, ti3))); |
| cr5 = VADD(SVMUL(ti11, tr5), SVMUL(ti12, tr4)); |
| ci5 = VADD(SVMUL(ti11, ti5), SVMUL(ti12, ti4)); |
| cr4 = VSUB(SVMUL(ti12, tr5), SVMUL(ti11, tr4)); |
| ci4 = VSUB(SVMUL(ti12, ti5), SVMUL(ti11, ti4)); |
| dr3 = VSUB(cr3, ci4); |
| dr4 = VADD(cr3, ci4); |
| di3 = VADD(ci3, cr4); |
| di4 = VSUB(ci3, cr4); |
| dr5 = VADD(cr2, ci5); |
| dr2 = VSUB(cr2, ci5); |
| di5 = VSUB(ci2, cr5); |
| di2 = VADD(ci2, cr5); |
| wr1=wa1[i], wi1=fsign*wa1[i+1], wr2=wa2[i], wi2=fsign*wa2[i+1]; |
| wr3=wa3[i], wi3=fsign*wa3[i+1], wr4=wa4[i], wi4=fsign*wa4[i+1]; |
| VCPLXMUL(dr2, di2, LD_PS1(wr1), LD_PS1(wi1)); |
| ch_ref(i - 1, 2) = dr2; |
| ch_ref(i, 2) = di2; |
| VCPLXMUL(dr3, di3, LD_PS1(wr2), LD_PS1(wi2)); |
| ch_ref(i - 1, 3) = dr3; |
| ch_ref(i, 3) = di3; |
| VCPLXMUL(dr4, di4, LD_PS1(wr3), LD_PS1(wi3)); |
| ch_ref(i - 1, 4) = dr4; |
| ch_ref(i, 4) = di4; |
| VCPLXMUL(dr5, di5, LD_PS1(wr4), LD_PS1(wi4)); |
| ch_ref(i - 1, 5) = dr5; |
| ch_ref(i, 5) = di5; |
| } |
| } |
| #undef ch_ref |
| #undef cc_ref |
| } |
| |
| static NEVER_INLINE(void) radf2_ps(int ido, int l1, const v4sf * RESTRICT cc, v4sf * RESTRICT ch, const float *wa1) { |
| static const float minus_one = -1.f; |
| int i, k, l1ido = l1*ido; |
| for (k=0; k < l1ido; k += ido) { |
| v4sf a = cc[k], b = cc[k + l1ido]; |
| ch[2*k] = VADD(a, b); |
| ch[2*(k+ido)-1] = VSUB(a, b); |
| } |
| if (ido < 2) return; |
| if (ido != 2) { |
| for (k=0; k < l1ido; k += ido) { |
| for (i=2; i<ido; i+=2) { |
| v4sf tr2 = cc[i - 1 + k + l1ido], ti2 = cc[i + k + l1ido]; |
| v4sf br = cc[i - 1 + k], bi = cc[i + k]; |
| VCPLXMULCONJ(tr2, ti2, LD_PS1(wa1[i - 2]), LD_PS1(wa1[i - 1])); |
| ch[i + 2*k] = VADD(bi, ti2); |
| ch[2*(k+ido) - i] = VSUB(ti2, bi); |
| ch[i - 1 + 2*k] = VADD(br, tr2); |
| ch[2*(k+ido) - i -1] = VSUB(br, tr2); |
| } |
| } |
| if (ido % 2 == 1) return; |
| } |
| for (k=0; k < l1ido; k += ido) { |
| ch[2*k + ido] = SVMUL(minus_one, cc[ido-1 + k + l1ido]); |
| ch[2*k + ido-1] = cc[k + ido-1]; |
| } |
| } /* radf2 */ |
| |
| |
| static NEVER_INLINE(void) radb2_ps(int ido, int l1, const v4sf *cc, v4sf *ch, const float *wa1) { |
| static const float minus_two=-2; |
| int i, k, l1ido = l1*ido; |
| v4sf a,b,c,d, tr2, ti2; |
| for (k=0; k < l1ido; k += ido) { |
| a = cc[2*k]; b = cc[2*(k+ido) - 1]; |
| ch[k] = VADD(a, b); |
| ch[k + l1ido] =VSUB(a, b); |
| } |
| if (ido < 2) return; |
| if (ido != 2) { |
| for (k = 0; k < l1ido; k += ido) { |
| for (i = 2; i < ido; i += 2) { |
| a = cc[i-1 + 2*k]; b = cc[2*(k + ido) - i - 1]; |
| c = cc[i+0 + 2*k]; d = cc[2*(k + ido) - i + 0]; |
| ch[i-1 + k] = VADD(a, b); |
| tr2 = VSUB(a, b); |
| ch[i+0 + k] = VSUB(c, d); |
| ti2 = VADD(c, d); |
| VCPLXMUL(tr2, ti2, LD_PS1(wa1[i - 2]), LD_PS1(wa1[i - 1])); |
| ch[i-1 + k + l1ido] = tr2; |
| ch[i+0 + k + l1ido] = ti2; |
| } |
| } |
| if (ido % 2 == 1) return; |
| } |
| for (k = 0; k < l1ido; k += ido) { |
| a = cc[2*k + ido-1]; b = cc[2*k + ido]; |
| ch[k + ido-1] = VADD(a,a); |
| ch[k + ido-1 + l1ido] = SVMUL(minus_two, b); |
| } |
| } /* radb2 */ |
| |
| static void radf3_ps(int ido, int l1, const v4sf * RESTRICT cc, v4sf * RESTRICT ch, |
| const float *wa1, const float *wa2) { |
| static const float taur = -0.5f; |
| static const float taui = 0.866025403784439f; |
| int i, k, ic; |
| v4sf ci2, di2, di3, cr2, dr2, dr3, ti2, ti3, tr2, tr3, wr1, wi1, wr2, wi2; |
| for (k=0; k<l1; k++) { |
| cr2 = VADD(cc[(k + l1)*ido], cc[(k + 2*l1)*ido]); |
| ch[3*k*ido] = VADD(cc[k*ido], cr2); |
| ch[(3*k+2)*ido] = SVMUL(taui, VSUB(cc[(k + l1*2)*ido], cc[(k + l1)*ido])); |
| ch[ido-1 + (3*k + 1)*ido] = VADD(cc[k*ido], SVMUL(taur, cr2)); |
| } |
| if (ido == 1) return; |
| for (k=0; k<l1; k++) { |
| for (i=2; i<ido; i+=2) { |
| ic = ido - i; |
| wr1 = LD_PS1(wa1[i - 2]); wi1 = LD_PS1(wa1[i - 1]); |
| dr2 = cc[i - 1 + (k + l1)*ido]; di2 = cc[i + (k + l1)*ido]; |
| VCPLXMULCONJ(dr2, di2, wr1, wi1); |
| |
| wr2 = LD_PS1(wa2[i - 2]); wi2 = LD_PS1(wa2[i - 1]); |
| dr3 = cc[i - 1 + (k + l1*2)*ido]; di3 = cc[i + (k + l1*2)*ido]; |
| VCPLXMULCONJ(dr3, di3, wr2, wi2); |
| |
| cr2 = VADD(dr2, dr3); |
| ci2 = VADD(di2, di3); |
| ch[i - 1 + 3*k*ido] = VADD(cc[i - 1 + k*ido], cr2); |
| ch[i + 3*k*ido] = VADD(cc[i + k*ido], ci2); |
| tr2 = VADD(cc[i - 1 + k*ido], SVMUL(taur, cr2)); |
| ti2 = VADD(cc[i + k*ido], SVMUL(taur, ci2)); |
| tr3 = SVMUL(taui, VSUB(di2, di3)); |
| ti3 = SVMUL(taui, VSUB(dr3, dr2)); |
| ch[i - 1 + (3*k + 2)*ido] = VADD(tr2, tr3); |
| ch[ic - 1 + (3*k + 1)*ido] = VSUB(tr2, tr3); |
| ch[i + (3*k + 2)*ido] = VADD(ti2, ti3); |
| ch[ic + (3*k + 1)*ido] = VSUB(ti3, ti2); |
| } |
| } |
| } /* radf3 */ |
| |
| |
| static void radb3_ps(int ido, int l1, const v4sf *RESTRICT cc, v4sf *RESTRICT ch, |
| const float *wa1, const float *wa2) |
| { |
| static const float taur = -0.5f; |
| static const float taui = 0.866025403784439f; |
| static const float taui_2 = 0.866025403784439f*2; |
| int i, k, ic; |
| v4sf ci2, ci3, di2, di3, cr2, cr3, dr2, dr3, ti2, tr2; |
| for (k=0; k<l1; k++) { |
| tr2 = cc[ido-1 + (3*k + 1)*ido]; tr2 = VADD(tr2,tr2); |
| cr2 = VMADD(LD_PS1(taur), tr2, cc[3*k*ido]); |
| ch[k*ido] = VADD(cc[3*k*ido], tr2); |
| ci3 = SVMUL(taui_2, cc[(3*k + 2)*ido]); |
| ch[(k + l1)*ido] = VSUB(cr2, ci3); |
| ch[(k + 2*l1)*ido] = VADD(cr2, ci3); |
| } |
| if (ido == 1) return; |
| for (k=0; k<l1; k++) { |
| for (i=2; i<ido; i+=2) { |
| ic = ido - i; |
| tr2 = VADD(cc[i - 1 + (3*k + 2)*ido], cc[ic - 1 + (3*k + 1)*ido]); |
| cr2 = VMADD(LD_PS1(taur), tr2, cc[i - 1 + 3*k*ido]); |
| ch[i - 1 + k*ido] = VADD(cc[i - 1 + 3*k*ido], tr2); |
| ti2 = VSUB(cc[i + (3*k + 2)*ido], cc[ic + (3*k + 1)*ido]); |
| ci2 = VMADD(LD_PS1(taur), ti2, cc[i + 3*k*ido]); |
| ch[i + k*ido] = VADD(cc[i + 3*k*ido], ti2); |
| cr3 = SVMUL(taui, VSUB(cc[i - 1 + (3*k + 2)*ido], cc[ic - 1 + (3*k + 1)*ido])); |
| ci3 = SVMUL(taui, VADD(cc[i + (3*k + 2)*ido], cc[ic + (3*k + 1)*ido])); |
| dr2 = VSUB(cr2, ci3); |
| dr3 = VADD(cr2, ci3); |
| di2 = VADD(ci2, cr3); |
| di3 = VSUB(ci2, cr3); |
| VCPLXMUL(dr2, di2, LD_PS1(wa1[i-2]), LD_PS1(wa1[i-1])); |
| ch[i - 1 + (k + l1)*ido] = dr2; |
| ch[i + (k + l1)*ido] = di2; |
| VCPLXMUL(dr3, di3, LD_PS1(wa2[i-2]), LD_PS1(wa2[i-1])); |
| ch[i - 1 + (k + 2*l1)*ido] = dr3; |
| ch[i + (k + 2*l1)*ido] = di3; |
| } |
| } |
| } /* radb3 */ |
| |
| static NEVER_INLINE(void) radf4_ps(int ido, int l1, const v4sf *RESTRICT cc, v4sf * RESTRICT ch, |
| const float * RESTRICT wa1, const float * RESTRICT wa2, const float * RESTRICT wa3) |
| { |
| static const float minus_hsqt2 = (float)-0.7071067811865475; |
| int i, k, l1ido = l1*ido; |
| { |
| const v4sf *RESTRICT cc_ = cc, * RESTRICT cc_end = cc + l1ido; |
| v4sf * RESTRICT ch_ = ch; |
| while (cc < cc_end) { |
| // this loop represents between 25% and 40% of total radf4_ps cost ! |
| v4sf a0 = cc[0], a1 = cc[l1ido]; |
| v4sf a2 = cc[2*l1ido], a3 = cc[3*l1ido]; |
| v4sf tr1 = VADD(a1, a3); |
| v4sf tr2 = VADD(a0, a2); |
| ch[2*ido-1] = VSUB(a0, a2); |
| ch[2*ido ] = VSUB(a3, a1); |
| ch[0 ] = VADD(tr1, tr2); |
| ch[4*ido-1] = VSUB(tr2, tr1); |
| cc += ido; ch += 4*ido; |
| } |
| cc = cc_; ch = ch_; |
| } |
| if (ido < 2) return; |
| if (ido != 2) { |
| for (k = 0; k < l1ido; k += ido) { |
| const v4sf * RESTRICT pc = (v4sf*)(cc + 1 + k); |
| for (i=2; i<ido; i += 2, pc += 2) { |
| int ic = ido - i; |
| v4sf wr, wi, cr2, ci2, cr3, ci3, cr4, ci4; |
| v4sf tr1, ti1, tr2, ti2, tr3, ti3, tr4, ti4; |
| |
| cr2 = pc[1*l1ido+0]; |
| ci2 = pc[1*l1ido+1]; |
| wr=LD_PS1(wa1[i - 2]); |
| wi=LD_PS1(wa1[i - 1]); |
| VCPLXMULCONJ(cr2,ci2,wr,wi); |
| |
| cr3 = pc[2*l1ido+0]; |
| ci3 = pc[2*l1ido+1]; |
| wr = LD_PS1(wa2[i-2]); |
| wi = LD_PS1(wa2[i-1]); |
| VCPLXMULCONJ(cr3, ci3, wr, wi); |
| |
| cr4 = pc[3*l1ido]; |
| ci4 = pc[3*l1ido+1]; |
| wr = LD_PS1(wa3[i-2]); |
| wi = LD_PS1(wa3[i-1]); |
| VCPLXMULCONJ(cr4, ci4, wr, wi); |
| |
| /* at this point, on SSE, five of "cr2 cr3 cr4 ci2 ci3 ci4" should be loaded in registers */ |
| |
| tr1 = VADD(cr2,cr4); |
| tr4 = VSUB(cr4,cr2); |
| tr2 = VADD(pc[0],cr3); |
| tr3 = VSUB(pc[0],cr3); |
| ch[i - 1 + 4*k] = VADD(tr1,tr2); |
| ch[ic - 1 + 4*k + 3*ido] = VSUB(tr2,tr1); // at this point tr1 and tr2 can be disposed |
| ti1 = VADD(ci2,ci4); |
| ti4 = VSUB(ci2,ci4); |
| ch[i - 1 + 4*k + 2*ido] = VADD(ti4,tr3); |
| ch[ic - 1 + 4*k + 1*ido] = VSUB(tr3,ti4); // dispose tr3, ti4 |
| ti2 = VADD(pc[1],ci3); |
| ti3 = VSUB(pc[1],ci3); |
| ch[i + 4*k] = VADD(ti1, ti2); |
| ch[ic + 4*k + 3*ido] = VSUB(ti1, ti2); |
| ch[i + 4*k + 2*ido] = VADD(tr4, ti3); |
| ch[ic + 4*k + 1*ido] = VSUB(tr4, ti3); |
| } |
| } |
| if (ido % 2 == 1) return; |
| } |
| for (k=0; k<l1ido; k += ido) { |
| v4sf a = cc[ido-1 + k + l1ido], b = cc[ido-1 + k + 3*l1ido]; |
| v4sf c = cc[ido-1 + k], d = cc[ido-1 + k + 2*l1ido]; |
| v4sf ti1 = SVMUL(minus_hsqt2, VADD(a, b)); |
| v4sf tr1 = SVMUL(minus_hsqt2, VSUB(b, a)); |
| ch[ido-1 + 4*k] = VADD(tr1, c); |
| ch[ido-1 + 4*k + 2*ido] = VSUB(c, tr1); |
| ch[4*k + 1*ido] = VSUB(ti1, d); |
| ch[4*k + 3*ido] = VADD(ti1, d); |
| } |
| } /* radf4 */ |
| |
| |
| static NEVER_INLINE(void) radb4_ps(int ido, int l1, const v4sf * RESTRICT cc, v4sf * RESTRICT ch, |
| const float * RESTRICT wa1, const float * RESTRICT wa2, const float *RESTRICT wa3) |
| { |
| static const float minus_sqrt2 = (float)-1.414213562373095; |
| static const float two = 2.f; |
| int i, k, l1ido = l1*ido; |
| v4sf ci2, ci3, ci4, cr2, cr3, cr4, ti1, ti2, ti3, ti4, tr1, tr2, tr3, tr4; |
| { |
| const v4sf *RESTRICT cc_ = cc, * RESTRICT ch_end = ch + l1ido; |
| v4sf *ch_ = ch; |
| while (ch < ch_end) { |
| v4sf a = cc[0], b = cc[4*ido-1]; |
| v4sf c = cc[2*ido], d = cc[2*ido-1]; |
| tr3 = SVMUL(two,d); |
| tr2 = VADD(a,b); |
| tr1 = VSUB(a,b); |
| tr4 = SVMUL(two,c); |
| ch[0*l1ido] = VADD(tr2, tr3); |
| ch[2*l1ido] = VSUB(tr2, tr3); |
| ch[1*l1ido] = VSUB(tr1, tr4); |
| ch[3*l1ido] = VADD(tr1, tr4); |
| |
| cc += 4*ido; ch += ido; |
| } |
| cc = cc_; ch = ch_; |
| } |
| if (ido < 2) return; |
| if (ido != 2) { |
| for (k = 0; k < l1ido; k += ido) { |
| const v4sf * RESTRICT pc = (v4sf*)(cc - 1 + 4*k); |
| v4sf * RESTRICT ph = (v4sf*)(ch + k + 1); |
| for (i = 2; i < ido; i += 2) { |
| |
| tr1 = VSUB(pc[i], pc[4*ido - i]); |
| tr2 = VADD(pc[i], pc[4*ido - i]); |
| ti4 = VSUB(pc[2*ido + i], pc[2*ido - i]); |
| tr3 = VADD(pc[2*ido + i], pc[2*ido - i]); |
| ph[0] = VADD(tr2, tr3); |
| cr3 = VSUB(tr2, tr3); |
| |
| ti3 = VSUB(pc[2*ido + i + 1], pc[2*ido - i + 1]); |
| tr4 = VADD(pc[2*ido + i + 1], pc[2*ido - i + 1]); |
| cr2 = VSUB(tr1, tr4); |
| cr4 = VADD(tr1, tr4); |
| |
| ti1 = VADD(pc[i + 1], pc[4*ido - i + 1]); |
| ti2 = VSUB(pc[i + 1], pc[4*ido - i + 1]); |
| |
| ph[1] = VADD(ti2, ti3); ph += l1ido; |
| ci3 = VSUB(ti2, ti3); |
| ci2 = VADD(ti1, ti4); |
| ci4 = VSUB(ti1, ti4); |
| VCPLXMUL(cr2, ci2, LD_PS1(wa1[i-2]), LD_PS1(wa1[i-1])); |
| ph[0] = cr2; |
| ph[1] = ci2; ph += l1ido; |
| VCPLXMUL(cr3, ci3, LD_PS1(wa2[i-2]), LD_PS1(wa2[i-1])); |
| ph[0] = cr3; |
| ph[1] = ci3; ph += l1ido; |
| VCPLXMUL(cr4, ci4, LD_PS1(wa3[i-2]), LD_PS1(wa3[i-1])); |
| ph[0] = cr4; |
| ph[1] = ci4; ph = ph - 3*l1ido + 2; |
| } |
| } |
| if (ido % 2 == 1) return; |
| } |
| for (k=0; k < l1ido; k+=ido) { |
| int i0 = 4*k + ido; |
| v4sf c = cc[i0-1], d = cc[i0 + 2*ido-1]; |
| v4sf a = cc[i0+0], b = cc[i0 + 2*ido+0]; |
| tr1 = VSUB(c,d); |
| tr2 = VADD(c,d); |
| ti1 = VADD(b,a); |
| ti2 = VSUB(b,a); |
| ch[ido-1 + k + 0*l1ido] = VADD(tr2,tr2); |
| ch[ido-1 + k + 1*l1ido] = SVMUL(minus_sqrt2, VSUB(ti1, tr1)); |
| ch[ido-1 + k + 2*l1ido] = VADD(ti2, ti2); |
| ch[ido-1 + k + 3*l1ido] = SVMUL(minus_sqrt2, VADD(ti1, tr1)); |
| } |
| } /* radb4 */ |
| |
| static void radf5_ps(int ido, int l1, const v4sf * RESTRICT cc, v4sf * RESTRICT ch, |
| const float *wa1, const float *wa2, const float *wa3, const float *wa4) |
| { |
| static const float tr11 = .309016994374947f; |
| static const float ti11 = .951056516295154f; |
| static const float tr12 = -.809016994374947f; |
| static const float ti12 = .587785252292473f; |
| |
| /* System generated locals */ |
| int cc_offset, ch_offset; |
| |
| /* Local variables */ |
| int i, k, ic; |
| v4sf ci2, di2, ci4, ci5, di3, di4, di5, ci3, cr2, cr3, dr2, dr3, dr4, dr5, |
| cr5, cr4, ti2, ti3, ti5, ti4, tr2, tr3, tr4, tr5; |
| int idp2; |
| |
| |
| #define cc_ref(a_1,a_2,a_3) cc[((a_3)*l1 + (a_2))*ido + a_1] |
| #define ch_ref(a_1,a_2,a_3) ch[((a_3)*5 + (a_2))*ido + a_1] |
| |
| /* Parameter adjustments */ |
| ch_offset = 1 + ido * 6; |
| ch -= ch_offset; |
| cc_offset = 1 + ido * (1 + l1); |
| cc -= cc_offset; |
| |
| /* Function Body */ |
| for (k = 1; k <= l1; ++k) { |
| cr2 = VADD(cc_ref(1, k, 5), cc_ref(1, k, 2)); |
| ci5 = VSUB(cc_ref(1, k, 5), cc_ref(1, k, 2)); |
| cr3 = VADD(cc_ref(1, k, 4), cc_ref(1, k, 3)); |
| ci4 = VSUB(cc_ref(1, k, 4), cc_ref(1, k, 3)); |
| ch_ref(1, 1, k) = VADD(cc_ref(1, k, 1), VADD(cr2, cr3)); |
| ch_ref(ido, 2, k) = VADD(cc_ref(1, k, 1), VADD(SVMUL(tr11, cr2), SVMUL(tr12, cr3))); |
| ch_ref(1, 3, k) = VADD(SVMUL(ti11, ci5), SVMUL(ti12, ci4)); |
| ch_ref(ido, 4, k) = VADD(cc_ref(1, k, 1), VADD(SVMUL(tr12, cr2), SVMUL(tr11, cr3))); |
| ch_ref(1, 5, k) = VSUB(SVMUL(ti12, ci5), SVMUL(ti11, ci4)); |
| //printf("pffft: radf5, k=%d ch_ref=%f, ci4=%f\n", k, ch_ref(1, 5, k), ci4); |
| } |
| if (ido == 1) { |
| return; |
| } |
| idp2 = ido + 2; |
| for (k = 1; k <= l1; ++k) { |
| for (i = 3; i <= ido; i += 2) { |
| ic = idp2 - i; |
| dr2 = LD_PS1(wa1[i-3]); di2 = LD_PS1(wa1[i-2]); |
| dr3 = LD_PS1(wa2[i-3]); di3 = LD_PS1(wa2[i-2]); |
| dr4 = LD_PS1(wa3[i-3]); di4 = LD_PS1(wa3[i-2]); |
| dr5 = LD_PS1(wa4[i-3]); di5 = LD_PS1(wa4[i-2]); |
| VCPLXMULCONJ(dr2, di2, cc_ref(i-1, k, 2), cc_ref(i, k, 2)); |
| VCPLXMULCONJ(dr3, di3, cc_ref(i-1, k, 3), cc_ref(i, k, 3)); |
| VCPLXMULCONJ(dr4, di4, cc_ref(i-1, k, 4), cc_ref(i, k, 4)); |
| VCPLXMULCONJ(dr5, di5, cc_ref(i-1, k, 5), cc_ref(i, k, 5)); |
| cr2 = VADD(dr2, dr5); |
| ci5 = VSUB(dr5, dr2); |
| cr5 = VSUB(di2, di5); |
| ci2 = VADD(di2, di5); |
| cr3 = VADD(dr3, dr4); |
| ci4 = VSUB(dr4, dr3); |
| cr4 = VSUB(di3, di4); |
| ci3 = VADD(di3, di4); |
| ch_ref(i - 1, 1, k) = VADD(cc_ref(i - 1, k, 1), VADD(cr2, cr3)); |
| ch_ref(i, 1, k) = VSUB(cc_ref(i, k, 1), VADD(ci2, ci3));// |
| tr2 = VADD(cc_ref(i - 1, k, 1), VADD(SVMUL(tr11, cr2), SVMUL(tr12, cr3))); |
| ti2 = VSUB(cc_ref(i, k, 1), VADD(SVMUL(tr11, ci2), SVMUL(tr12, ci3)));// |
| tr3 = VADD(cc_ref(i - 1, k, 1), VADD(SVMUL(tr12, cr2), SVMUL(tr11, cr3))); |
| ti3 = VSUB(cc_ref(i, k, 1), VADD(SVMUL(tr12, ci2), SVMUL(tr11, ci3)));// |
| tr5 = VADD(SVMUL(ti11, cr5), SVMUL(ti12, cr4)); |
| ti5 = VADD(SVMUL(ti11, ci5), SVMUL(ti12, ci4)); |
| tr4 = VSUB(SVMUL(ti12, cr5), SVMUL(ti11, cr4)); |
| ti4 = VSUB(SVMUL(ti12, ci5), SVMUL(ti11, ci4)); |
| ch_ref(i - 1, 3, k) = VSUB(tr2, tr5); |
| ch_ref(ic - 1, 2, k) = VADD(tr2, tr5); |
| ch_ref(i, 3, k) = VADD(ti2, ti5); |
| ch_ref(ic, 2, k) = VSUB(ti5, ti2); |
| ch_ref(i - 1, 5, k) = VSUB(tr3, tr4); |
| ch_ref(ic - 1, 4, k) = VADD(tr3, tr4); |
| ch_ref(i, 5, k) = VADD(ti3, ti4); |
| ch_ref(ic, 4, k) = VSUB(ti4, ti3); |
| } |
| } |
| #undef cc_ref |
| #undef ch_ref |
| } /* radf5 */ |
| |
| static void radb5_ps(int ido, int l1, const v4sf *RESTRICT cc, v4sf *RESTRICT ch, |
| const float *wa1, const float *wa2, const float *wa3, const float *wa4) |
| { |
| static const float tr11 = .309016994374947f; |
| static const float ti11 = .951056516295154f; |
| static const float tr12 = -.809016994374947f; |
| static const float ti12 = .587785252292473f; |
| |
| int cc_offset, ch_offset; |
| |
| /* Local variables */ |
| int i, k, ic; |
| v4sf ci2, ci3, ci4, ci5, di3, di4, di5, di2, cr2, cr3, cr5, cr4, ti2, ti3, |
| ti4, ti5, dr3, dr4, dr5, dr2, tr2, tr3, tr4, tr5; |
| int idp2; |
| |
| #define cc_ref(a_1,a_2,a_3) cc[((a_3)*5 + (a_2))*ido + a_1] |
| #define ch_ref(a_1,a_2,a_3) ch[((a_3)*l1 + (a_2))*ido + a_1] |
| |
| /* Parameter adjustments */ |
| ch_offset = 1 + ido * (1 + l1); |
| ch -= ch_offset; |
| cc_offset = 1 + ido * 6; |
| cc -= cc_offset; |
| |
| /* Function Body */ |
| for (k = 1; k <= l1; ++k) { |
| ti5 = VADD(cc_ref(1, 3, k), cc_ref(1, 3, k)); |
| ti4 = VADD(cc_ref(1, 5, k), cc_ref(1, 5, k)); |
| tr2 = VADD(cc_ref(ido, 2, k), cc_ref(ido, 2, k)); |
| tr3 = VADD(cc_ref(ido, 4, k), cc_ref(ido, 4, k)); |
| ch_ref(1, k, 1) = VADD(cc_ref(1, 1, k), VADD(tr2, tr3)); |
| cr2 = VADD(cc_ref(1, 1, k), VADD(SVMUL(tr11, tr2), SVMUL(tr12, tr3))); |
| cr3 = VADD(cc_ref(1, 1, k), VADD(SVMUL(tr12, tr2), SVMUL(tr11, tr3))); |
| ci5 = VADD(SVMUL(ti11, ti5), SVMUL(ti12, ti4)); |
| ci4 = VSUB(SVMUL(ti12, ti5), SVMUL(ti11, ti4)); |
| ch_ref(1, k, 2) = VSUB(cr2, ci5); |
| ch_ref(1, k, 3) = VSUB(cr3, ci4); |
| ch_ref(1, k, 4) = VADD(cr3, ci4); |
| ch_ref(1, k, 5) = VADD(cr2, ci5); |
| } |
| if (ido == 1) { |
| return; |
| } |
| idp2 = ido + 2; |
| for (k = 1; k <= l1; ++k) { |
| for (i = 3; i <= ido; i += 2) { |
| ic = idp2 - i; |
| ti5 = VADD(cc_ref(i , 3, k), cc_ref(ic , 2, k)); |
| ti2 = VSUB(cc_ref(i , 3, k), cc_ref(ic , 2, k)); |
| ti4 = VADD(cc_ref(i , 5, k), cc_ref(ic , 4, k)); |
| ti3 = VSUB(cc_ref(i , 5, k), cc_ref(ic , 4, k)); |
| tr5 = VSUB(cc_ref(i-1, 3, k), cc_ref(ic-1, 2, k)); |
| tr2 = VADD(cc_ref(i-1, 3, k), cc_ref(ic-1, 2, k)); |
| tr4 = VSUB(cc_ref(i-1, 5, k), cc_ref(ic-1, 4, k)); |
| tr3 = VADD(cc_ref(i-1, 5, k), cc_ref(ic-1, 4, k)); |
| ch_ref(i - 1, k, 1) = VADD(cc_ref(i-1, 1, k), VADD(tr2, tr3)); |
| ch_ref(i, k, 1) = VADD(cc_ref(i, 1, k), VADD(ti2, ti3)); |
| cr2 = VADD(cc_ref(i-1, 1, k), VADD(SVMUL(tr11, tr2), SVMUL(tr12, tr3))); |
| ci2 = VADD(cc_ref(i , 1, k), VADD(SVMUL(tr11, ti2), SVMUL(tr12, ti3))); |
| cr3 = VADD(cc_ref(i-1, 1, k), VADD(SVMUL(tr12, tr2), SVMUL(tr11, tr3))); |
| ci3 = VADD(cc_ref(i , 1, k), VADD(SVMUL(tr12, ti2), SVMUL(tr11, ti3))); |
| cr5 = VADD(SVMUL(ti11, tr5), SVMUL(ti12, tr4)); |
| ci5 = VADD(SVMUL(ti11, ti5), SVMUL(ti12, ti4)); |
| cr4 = VSUB(SVMUL(ti12, tr5), SVMUL(ti11, tr4)); |
| ci4 = VSUB(SVMUL(ti12, ti5), SVMUL(ti11, ti4)); |
| dr3 = VSUB(cr3, ci4); |
| dr4 = VADD(cr3, ci4); |
| di3 = VADD(ci3, cr4); |
| di4 = VSUB(ci3, cr4); |
| dr5 = VADD(cr2, ci5); |
| dr2 = VSUB(cr2, ci5); |
| di5 = VSUB(ci2, cr5); |
| di2 = VADD(ci2, cr5); |
| VCPLXMUL(dr2, di2, LD_PS1(wa1[i-3]), LD_PS1(wa1[i-2])); |
| VCPLXMUL(dr3, di3, LD_PS1(wa2[i-3]), LD_PS1(wa2[i-2])); |
| VCPLXMUL(dr4, di4, LD_PS1(wa3[i-3]), LD_PS1(wa3[i-2])); |
| VCPLXMUL(dr5, di5, LD_PS1(wa4[i-3]), LD_PS1(wa4[i-2])); |
| |
| ch_ref(i-1, k, 2) = dr2; ch_ref(i, k, 2) = di2; |
| ch_ref(i-1, k, 3) = dr3; ch_ref(i, k, 3) = di3; |
| ch_ref(i-1, k, 4) = dr4; ch_ref(i, k, 4) = di4; |
| ch_ref(i-1, k, 5) = dr5; ch_ref(i, k, 5) = di5; |
| } |
| } |
| #undef cc_ref |
| #undef ch_ref |
| } /* radb5 */ |
| |
| static NEVER_INLINE(v4sf *) rfftf1_ps(int n, const v4sf *input_readonly, v4sf *work1, v4sf *work2, |
| const float *wa, const int *ifac) { |
| v4sf *in = (v4sf*)input_readonly; |
| v4sf *out = (in == work2 ? work1 : work2); |
| int nf = ifac[1], k1; |
| int l2 = n; |
| int iw = n-1; |
| assert(in != out && work1 != work2); |
| for (k1 = 1; k1 <= nf; ++k1) { |
| int kh = nf - k1; |
| int ip = ifac[kh + 2]; |
| int l1 = l2 / ip; |
| int ido = n / l2; |
| iw -= (ip - 1)*ido; |
| switch (ip) { |
| case 5: { |
| int ix2 = iw + ido; |
| int ix3 = ix2 + ido; |
| int ix4 = ix3 + ido; |
| radf5_ps(ido, l1, in, out, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); |
| } break; |
| case 4: { |
| int ix2 = iw + ido; |
| int ix3 = ix2 + ido; |
| radf4_ps(ido, l1, in, out, &wa[iw], &wa[ix2], &wa[ix3]); |
| } break; |
| case 3: { |
| int ix2 = iw + ido; |
| radf3_ps(ido, l1, in, out, &wa[iw], &wa[ix2]); |
| } break; |
| case 2: |
| radf2_ps(ido, l1, in, out, &wa[iw]); |
| break; |
| default: |
| assert(0); |
| break; |
| } |
| l2 = l1; |
| if (out == work2) { |
| out = work1; in = work2; |
| } else { |
| out = work2; in = work1; |
| } |
| } |
| return in; /* this is in fact the output .. */ |
| } /* rfftf1 */ |
| |
| static NEVER_INLINE(v4sf *) rfftb1_ps(int n, const v4sf *input_readonly, v4sf *work1, v4sf *work2, |
| const float *wa, const int *ifac) { |
| v4sf *in = (v4sf*)input_readonly; |
| v4sf *out = (in == work2 ? work1 : work2); |
| int nf = ifac[1], k1; |
| int l1 = 1; |
| int iw = 0; |
| assert(in != out); |
| for (k1=1; k1<=nf; k1++) { |
| int ip = ifac[k1 + 1]; |
| int l2 = ip*l1; |
| int ido = n / l2; |
| switch (ip) { |
| case 5: { |
| int ix2 = iw + ido; |
| int ix3 = ix2 + ido; |
| int ix4 = ix3 + ido; |
| radb5_ps(ido, l1, in, out, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4]); |
| } break; |
| case 4: { |
| int ix2 = iw + ido; |
| int ix3 = ix2 + ido; |
| radb4_ps(ido, l1, in, out, &wa[iw], &wa[ix2], &wa[ix3]); |
| } break; |
| case 3: { |
| int ix2 = iw + ido; |
| radb3_ps(ido, l1, in, out, &wa[iw], &wa[ix2]); |
| } break; |
| case 2: |
| radb2_ps(ido, l1, in, out, &wa[iw]); |
| break; |
| default: |
| assert(0); |
| break; |
| } |
| l1 = l2; |
| iw += (ip - 1)*ido; |
| |
| if (out == work2) { |
| out = work1; in = work2; |
| } else { |
| out = work2; in = work1; |
| } |
| } |
| return in; /* this is in fact the output .. */ |
| } |
| |
| static int decompose(int n, int *ifac, const int *ntryh) { |
| int nl = n, nf = 0, i, j = 0; |
| for (j=0; ntryh[j]; ++j) { |
| int ntry = ntryh[j]; |
| while (nl != 1) { |
| int nq = nl / ntry; |
| int nr = nl - ntry * nq; |
| if (nr == 0) { |
| ifac[2+nf++] = ntry; |
| nl = nq; |
| if (ntry == 2 && nf != 1) { |
| for (i = 2; i <= nf; ++i) { |
| int ib = nf - i + 2; |
| ifac[ib + 1] = ifac[ib]; |
| } |
| ifac[2] = 2; |
| } |
| } else break; |
| } |
| } |
| ifac[0] = n; |
| ifac[1] = nf; |
| return nf; |
| } |
| |
| |
| |
| static void rffti1_ps(int n, float *wa, int *ifac) |
| { |
| static const int ntryh[] = { 4,2,3,5,0 }; |
| int k1, j, ii; |
| |
| int nf = decompose(n,ifac,ntryh); |
| float argh = (2*(float)M_PI) / n; |
| int is = 0; |
| int nfm1 = nf - 1; |
| int l1 = 1; |
| for (k1 = 1; k1 <= nfm1; k1++) { |
| int ip = ifac[k1 + 1]; |
| int ld = 0; |
| int l2 = l1*ip; |
| int ido = n / l2; |
| int ipm = ip - 1; |
| for (j = 1; j <= ipm; ++j) { |
| float argld; |
| int i = is, fi=0; |
| ld += l1; |
| argld = ld*argh; |
| for (ii = 3; ii <= ido; ii += 2) { |
| i += 2; |
| fi += 1; |
| wa[i - 2] = cosf(fi*argld); |
| wa[i - 1] = sinf(fi*argld); |
| } |
| is += ido; |
| } |
| l1 = l2; |
| } |
| } /* rffti1 */ |
| |
| void cffti1_ps(int n, float *wa, int *ifac) |
| { |
| static const int ntryh[] = { 5,3,4,2,0 }; |
| int k1, j, ii; |
| |
| int nf = decompose(n,ifac,ntryh); |
| float argh = (2*(float)M_PI) / n; |
| int i = 1; |
| int l1 = 1; |
| for (k1=1; k1<=nf; k1++) { |
| int ip = ifac[k1+1]; |
| int ld = 0; |
| int l2 = l1*ip; |
| int ido = n / l2; |
| int idot = ido + ido + 2; |
| int ipm = ip - 1; |
| for (j=1; j<=ipm; j++) { |
| float argld; |
| int i1 = i, fi = 0; |
| wa[i-1] = 1; |
| wa[i] = 0; |
| ld += l1; |
| argld = ld*argh; |
| for (ii = 4; ii <= idot; ii += 2) { |
| i += 2; |
| fi += 1; |
| wa[i-1] = cosf(fi*argld); |
| wa[i] = sinf(fi*argld); |
| } |
| if (ip > 5) { |
| wa[i1-1] = wa[i-1]; |
| wa[i1] = wa[i]; |
| } |
| } |
| l1 = l2; |
| } |
| } /* cffti1 */ |
| |
| |
| v4sf *cfftf1_ps(int n, const v4sf *input_readonly, v4sf *work1, v4sf *work2, const float *wa, const int *ifac, int isign) { |
| v4sf *in = (v4sf*)input_readonly; |
| v4sf *out = (in == work2 ? work1 : work2); |
| int nf = ifac[1], k1; |
| int l1 = 1; |
| int iw = 0; |
| assert(in != out && work1 != work2); |
| for (k1=2; k1<=nf+1; k1++) { |
| int ip = ifac[k1]; |
| int l2 = ip*l1; |
| int ido = n / l2; |
| int idot = ido + ido; |
| switch (ip) { |
| case 5: { |
| int ix2 = iw + idot; |
| int ix3 = ix2 + idot; |
| int ix4 = ix3 + idot; |
| passf5_ps(idot, l1, in, out, &wa[iw], &wa[ix2], &wa[ix3], &wa[ix4], isign); |
| } break; |
| case 4: { |
| int ix2 = iw + idot; |
| int ix3 = ix2 + idot; |
| passf4_ps(idot, l1, in, out, &wa[iw], &wa[ix2], &wa[ix3], isign); |
| } break; |
| case 2: { |
| passf2_ps(idot, l1, in, out, &wa[iw], isign); |
| } break; |
| case 3: { |
| int ix2 = iw + idot; |
| passf3_ps(idot, l1, in, out, &wa[iw], &wa[ix2], isign); |
| } break; |
| default: |
| assert(0); |
| } |
| l1 = l2; |
| iw += (ip - 1)*idot; |
| if (out == work2) { |
| out = work1; in = work2; |
| } else { |
| out = work2; in = work1; |
| } |
| } |
| |
| return in; /* this is in fact the output .. */ |
| } |
| |
| |
| struct PFFFT_Setup { |
| int N; |
| int Ncvec; // nb of complex simd vectors (N/4 if PFFFT_COMPLEX, N/8 if PFFFT_REAL) |
| int ifac[15]; |
| pffft_transform_t transform; |
| v4sf *data; // allocated room for twiddle coefs |
| float *e; // points into 'data' , N/4*3 elements |
| float *twiddle; // points into 'data', N/4 elements |
| }; |
| |
| PFFFT_Setup *pffft_new_setup(int N, pffft_transform_t transform) { |
| PFFFT_Setup *s = (PFFFT_Setup*)malloc(sizeof(PFFFT_Setup)); |
| int k, m; |
| /* unfortunately, the fft size must be a multiple of 16 for complex FFTs |
| and 32 for real FFTs -- a lot of stuff would need to be rewritten to |
| handle other cases (or maybe just switch to a scalar fft, I don't know..) */ |
| if (transform == PFFFT_REAL) { assert((N%(2*SIMD_SZ*SIMD_SZ))==0 && N>0); } |
| if (transform == PFFFT_COMPLEX) { assert((N%(SIMD_SZ*SIMD_SZ))==0 && N>0); } |
| //assert((N % 32) == 0); |
| s->N = N; |
| s->transform = transform; |
| /* nb of complex simd vectors */ |
| s->Ncvec = (transform == PFFFT_REAL ? N/2 : N)/SIMD_SZ; |
| s->data = (v4sf*)pffft_aligned_malloc(2*s->Ncvec * sizeof(v4sf)); |
| s->e = (float*)s->data; |
| s->twiddle = (float*)(s->data + (2*s->Ncvec*(SIMD_SZ-1))/SIMD_SZ); |
| |
| if (transform == PFFFT_REAL) { |
| for (k=0; k < s->Ncvec; ++k) { |
| int i = k/SIMD_SZ; |
| int j = k%SIMD_SZ; |
| for (m=0; m < SIMD_SZ-1; ++m) { |
| float A = -2*(float)M_PI*(m+1)*k / N; |
| s->e[(2*(i*3 + m) + 0) * SIMD_SZ + j] = cosf(A); |
| s->e[(2*(i*3 + m) + 1) * SIMD_SZ + j] = sinf(A); |
| } |
| } |
| rffti1_ps(N/SIMD_SZ, s->twiddle, s->ifac); |
| } else { |
| for (k=0; k < s->Ncvec; ++k) { |
| int i = k/SIMD_SZ; |
| int j = k%SIMD_SZ; |
| for (m=0; m < SIMD_SZ-1; ++m) { |
| float A = -2*(float)M_PI*(m+1)*k / N; |
| s->e[(2*(i*3 + m) + 0)*SIMD_SZ + j] = cosf(A); |
| s->e[(2*(i*3 + m) + 1)*SIMD_SZ + j] = sinf(A); |
| } |
| } |
| cffti1_ps(N/SIMD_SZ, s->twiddle, s->ifac); |
| } |
| |
| /* check that N is decomposable with allowed prime factors */ |
| for (k=0, m=1; k < s->ifac[1]; ++k) { m *= s->ifac[2+k]; } |
| if (m != N/SIMD_SZ) { |
| pffft_destroy_setup(s); s = 0; |
| } |
| |
| return s; |
| } |
| |
| |
| void pffft_destroy_setup(PFFFT_Setup *s) { |
| pffft_aligned_free(s->data); |
| free(s); |
| } |
| |
| #if !defined(PFFFT_SIMD_DISABLE) |
| |
| /* [0 0 1 2 3 4 5 6 7 8] -> [0 8 7 6 5 4 3 2 1] */ |
| static void reversed_copy(int N, const v4sf *in, int in_stride, v4sf *out) { |
| v4sf g0, g1; |
| int k; |
| INTERLEAVE2(in[0], in[1], g0, g1); in += in_stride; |
| |
| *--out = VSWAPHL(g0, g1); // [g0l, g0h], [g1l g1h] -> [g1l, g0h] |
| for (k=1; k < N; ++k) { |
| v4sf h0, h1; |
| INTERLEAVE2(in[0], in[1], h0, h1); in += in_stride; |
| *--out = VSWAPHL(g1, h0); |
| *--out = VSWAPHL(h0, h1); |
| g1 = h1; |
| } |
| *--out = VSWAPHL(g1, g0); |
| } |
| |
| static void unreversed_copy(int N, const v4sf *in, v4sf *out, int out_stride) { |
| v4sf g0, g1, h0, h1; |
| int k; |
| g0 = g1 = in[0]; ++in; |
| for (k=1; k < N; ++k) { |
| h0 = *in++; h1 = *in++; |
| g1 = VSWAPHL(g1, h0); |
| h0 = VSWAPHL(h0, h1); |
| UNINTERLEAVE2(h0, g1, out[0], out[1]); out += out_stride; |
| g1 = h1; |
| } |
| h0 = *in++; h1 = g0; |
| g1 = VSWAPHL(g1, h0); |
| h0 = VSWAPHL(h0, h1); |
| UNINTERLEAVE2(h0, g1, out[0], out[1]); |
| } |
| |
| void pffft_zreorder(PFFFT_Setup *setup, const float *in, float *out, pffft_direction_t direction) { |
| int k, N = setup->N, Ncvec = setup->Ncvec; |
| const v4sf *vin = (const v4sf*)in; |
| v4sf *vout = (v4sf*)out; |
| assert(in != out); |
| if (setup->transform == PFFFT_REAL) { |
| int k, dk = N/32; |
| if (direction == PFFFT_FORWARD) { |
| for (k=0; k < dk; ++k) { |
| INTERLEAVE2(vin[k*8 + 0], vin[k*8 + 1], vout[2*(0*dk + k) + 0], vout[2*(0*dk + k) + 1]); |
| INTERLEAVE2(vin[k*8 + 4], vin[k*8 + 5], vout[2*(2*dk + k) + 0], vout[2*(2*dk + k) + 1]); |
| } |
| reversed_copy(dk, vin+2, 8, (v4sf*)(out + N/2)); |
| reversed_copy(dk, vin+6, 8, (v4sf*)(out + N)); |
| } else { |
| for (k=0; k < dk; ++k) { |
| UNINTERLEAVE2(vin[2*(0*dk + k) + 0], vin[2*(0*dk + k) + 1], vout[k*8 + 0], vout[k*8 + 1]); |
| UNINTERLEAVE2(vin[2*(2*dk + k) + 0], vin[2*(2*dk + k) + 1], vout[k*8 + 4], vout[k*8 + 5]); |
| } |
| unreversed_copy(dk, (v4sf*)(in + N/4), (v4sf*)(out + N - 6*SIMD_SZ), -8); |
| unreversed_copy(dk, (v4sf*)(in + 3*N/4), (v4sf*)(out + N - 2*SIMD_SZ), -8); |
| } |
| } else { |
| if (direction == PFFFT_FORWARD) { |
| for (k=0; k < Ncvec; ++k) { |
| int kk = (k/4) + (k%4)*(Ncvec/4); |
| INTERLEAVE2(vin[k*2], vin[k*2+1], vout[kk*2], vout[kk*2+1]); |
| } |
| } else { |
| for (k=0; k < Ncvec; ++k) { |
| int kk = (k/4) + (k%4)*(Ncvec/4); |
| UNINTERLEAVE2(vin[kk*2], vin[kk*2+1], vout[k*2], vout[k*2+1]); |
| } |
| } |
| } |
| } |
| |
| void pffft_cplx_finalize(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e) { |
| int k, dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks |
| v4sf r0, i0, r1, i1, r2, i2, r3, i3; |
| v4sf sr0, dr0, sr1, dr1, si0, di0, si1, di1; |
| assert(in != out); |
| for (k=0; k < dk; ++k) { |
| r0 = in[8*k+0]; i0 = in[8*k+1]; |
| r1 = in[8*k+2]; i1 = in[8*k+3]; |
| r2 = in[8*k+4]; i2 = in[8*k+5]; |
| r3 = in[8*k+6]; i3 = in[8*k+7]; |
| VTRANSPOSE4(r0,r1,r2,r3); |
| VTRANSPOSE4(i0,i1,i2,i3); |
| VCPLXMUL(r1,i1,e[k*6+0],e[k*6+1]); |
| VCPLXMUL(r2,i2,e[k*6+2],e[k*6+3]); |
| VCPLXMUL(r3,i3,e[k*6+4],e[k*6+5]); |
| |
| sr0 = VADD(r0,r2); dr0 = VSUB(r0, r2); |
| sr1 = VADD(r1,r3); dr1 = VSUB(r1, r3); |
| si0 = VADD(i0,i2); di0 = VSUB(i0, i2); |
| si1 = VADD(i1,i3); di1 = VSUB(i1, i3); |
| |
| /* |
| transformation for each column is: |
| |
| [1 1 1 1 0 0 0 0] [r0] |
| [1 0 -1 0 0 -1 0 1] [r1] |
| [1 -1 1 -1 0 0 0 0] [r2] |
| [1 0 -1 0 0 1 0 -1] [r3] |
| [0 0 0 0 1 1 1 1] * [i0] |
| [0 1 0 -1 1 0 -1 0] [i1] |
| [0 0 0 0 1 -1 1 -1] [i2] |
| [0 -1 0 1 1 0 -1 0] [i3] |
| */ |
| |
| r0 = VADD(sr0, sr1); i0 = VADD(si0, si1); |
| r1 = VADD(dr0, di1); i1 = VSUB(di0, dr1); |
| r2 = VSUB(sr0, sr1); i2 = VSUB(si0, si1); |
| r3 = VSUB(dr0, di1); i3 = VADD(di0, dr1); |
| |
| *out++ = r0; *out++ = i0; *out++ = r1; *out++ = i1; |
| *out++ = r2; *out++ = i2; *out++ = r3; *out++ = i3; |
| } |
| } |
| |
| void pffft_cplx_preprocess(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e) { |
| int k, dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks |
| v4sf r0, i0, r1, i1, r2, i2, r3, i3; |
| v4sf sr0, dr0, sr1, dr1, si0, di0, si1, di1; |
| assert(in != out); |
| for (k=0; k < dk; ++k) { |
| r0 = in[8*k+0]; i0 = in[8*k+1]; |
| r1 = in[8*k+2]; i1 = in[8*k+3]; |
| r2 = in[8*k+4]; i2 = in[8*k+5]; |
| r3 = in[8*k+6]; i3 = in[8*k+7]; |
| |
| sr0 = VADD(r0,r2); dr0 = VSUB(r0, r2); |
| sr1 = VADD(r1,r3); dr1 = VSUB(r1, r3); |
| si0 = VADD(i0,i2); di0 = VSUB(i0, i2); |
| si1 = VADD(i1,i3); di1 = VSUB(i1, i3); |
| |
| r0 = VADD(sr0, sr1); i0 = VADD(si0, si1); |
| r1 = VSUB(dr0, di1); i1 = VADD(di0, dr1); |
| r2 = VSUB(sr0, sr1); i2 = VSUB(si0, si1); |
| r3 = VADD(dr0, di1); i3 = VSUB(di0, dr1); |
| |
| VCPLXMULCONJ(r1,i1,e[k*6+0],e[k*6+1]); |
| VCPLXMULCONJ(r2,i2,e[k*6+2],e[k*6+3]); |
| VCPLXMULCONJ(r3,i3,e[k*6+4],e[k*6+5]); |
| |
| VTRANSPOSE4(r0,r1,r2,r3); |
| VTRANSPOSE4(i0,i1,i2,i3); |
| |
| *out++ = r0; *out++ = i0; *out++ = r1; *out++ = i1; |
| *out++ = r2; *out++ = i2; *out++ = r3; *out++ = i3; |
| } |
| } |
| |
| |
| static ALWAYS_INLINE(void) pffft_real_finalize_4x4(const v4sf *in0, const v4sf *in1, const v4sf *in, |
| const v4sf *e, v4sf *out) { |
| v4sf r0, i0, r1, i1, r2, i2, r3, i3; |
| v4sf sr0, dr0, sr1, dr1, si0, di0, si1, di1; |
| r0 = *in0; i0 = *in1; |
| r1 = *in++; i1 = *in++; r2 = *in++; i2 = *in++; r3 = *in++; i3 = *in++; |
| VTRANSPOSE4(r0,r1,r2,r3); |
| VTRANSPOSE4(i0,i1,i2,i3); |
| |
| /* |
| transformation for each column is: |
| |
| [1 1 1 1 0 0 0 0] [r0] |
| [1 0 -1 0 0 -1 0 1] [r1] |
| [1 0 -1 0 0 1 0 -1] [r2] |
| [1 -1 1 -1 0 0 0 0] [r3] |
| [0 0 0 0 1 1 1 1] * [i0] |
| [0 -1 0 1 -1 0 1 0] [i1] |
| [0 -1 0 1 1 0 -1 0] [i2] |
| [0 0 0 0 -1 1 -1 1] [i3] |
| */ |
| |
| //cerr << "matrix initial, before e , REAL:\n 1: " << r0 << "\n 1: " << r1 << "\n 1: " << r2 << "\n 1: " << r3 << "\n"; |
| //cerr << "matrix initial, before e, IMAG :\n 1: " << i0 << "\n 1: " << i1 << "\n 1: " << i2 << "\n 1: " << i3 << "\n"; |
| |
| VCPLXMUL(r1,i1,e[0],e[1]); |
| VCPLXMUL(r2,i2,e[2],e[3]); |
| VCPLXMUL(r3,i3,e[4],e[5]); |
| |
| //cerr << "matrix initial, real part:\n 1: " << r0 << "\n 1: " << r1 << "\n 1: " << r2 << "\n 1: " << r3 << "\n"; |
| //cerr << "matrix initial, imag part:\n 1: " << i0 << "\n 1: " << i1 << "\n 1: " << i2 << "\n 1: " << i3 << "\n"; |
| |
| sr0 = VADD(r0,r2); dr0 = VSUB(r0,r2); |
| sr1 = VADD(r1,r3); dr1 = VSUB(r3,r1); |
| si0 = VADD(i0,i2); di0 = VSUB(i0,i2); |
| si1 = VADD(i1,i3); di1 = VSUB(i3,i1); |
| |
| r0 = VADD(sr0, sr1); |
| r3 = VSUB(sr0, sr1); |
| i0 = VADD(si0, si1); |
| i3 = VSUB(si1, si0); |
| r1 = VADD(dr0, di1); |
| r2 = VSUB(dr0, di1); |
| i1 = VSUB(dr1, di0); |
| i2 = VADD(dr1, di0); |
| |
| *out++ = r0; |
| *out++ = i0; |
| *out++ = r1; |
| *out++ = i1; |
| *out++ = r2; |
| *out++ = i2; |
| *out++ = r3; |
| *out++ = i3; |
| |
| } |
| |
| static NEVER_INLINE(void) pffft_real_finalize(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e) { |
| int k, dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks |
| /* fftpack order is f0r f1r f1i f2r f2i ... f(n-1)r f(n-1)i f(n)r */ |
| |
| v4sf_union cr, ci, *uout = (v4sf_union*)out; |
| v4sf save = in[7], zero=VZERO(); |
| float xr0, xi0, xr1, xi1, xr2, xi2, xr3, xi3; |
| static const float s = (float)M_SQRT2/2; |
| |
| cr.v = in[0]; ci.v = in[Ncvec*2-1]; |
| assert(in != out); |
| pffft_real_finalize_4x4(&zero, &zero, in+1, e, out); |
| |
| /* |
| [cr0 cr1 cr2 cr3 ci0 ci1 ci2 ci3] |
| |
| [Xr(1)] ] [1 1 1 1 0 0 0 0] |
| [Xr(N/4) ] [0 0 0 0 1 s 0 -s] |
| [Xr(N/2) ] [1 0 -1 0 0 0 0 0] |
| [Xr(3N/4)] [0 0 0 0 1 -s 0 s] |
| [Xi(1) ] [1 -1 1 -1 0 0 0 0] |
| [Xi(N/4) ] [0 0 0 0 0 -s -1 -s] |
| [Xi(N/2) ] [0 -1 0 1 0 0 0 0] |
| [Xi(3N/4)] [0 0 0 0 0 -s 1 -s] |
| */ |
| |
| xr0=(cr.f[0]+cr.f[2]) + (cr.f[1]+cr.f[3]); uout[0].f[0] = xr0; |
| xi0=(cr.f[0]+cr.f[2]) - (cr.f[1]+cr.f[3]); uout[1].f[0] = xi0; |
| xr2=(cr.f[0]-cr.f[2]); uout[4].f[0] = xr2; |
| xi2=(cr.f[3]-cr.f[1]); uout[5].f[0] = xi2; |
| xr1= ci.f[0] + s*(ci.f[1]-ci.f[3]); uout[2].f[0] = xr1; |
| xi1=-ci.f[2] - s*(ci.f[1]+ci.f[3]); uout[3].f[0] = xi1; |
| xr3= ci.f[0] - s*(ci.f[1]-ci.f[3]); uout[6].f[0] = xr3; |
| xi3= ci.f[2] - s*(ci.f[1]+ci.f[3]); uout[7].f[0] = xi3; |
| |
| for (k=1; k < dk; ++k) { |
| v4sf save_next = in[8*k+7]; |
| pffft_real_finalize_4x4(&save, &in[8*k+0], in + 8*k+1, |
| e + k*6, out + k*8); |
| save = save_next; |
| } |
| |
| } |
| |
| static ALWAYS_INLINE(void) pffft_real_preprocess_4x4(const v4sf *in, |
| const v4sf *e, v4sf *out, int first) { |
| v4sf r0=in[0], i0=in[1], r1=in[2], i1=in[3], r2=in[4], i2=in[5], r3=in[6], i3=in[7]; |
| /* |
| transformation for each column is: |
| |
| [1 1 1 1 0 0 0 0] [r0] |
| [1 0 0 -1 0 -1 -1 0] [r1] |
| [1 -1 -1 1 0 0 0 0] [r2] |
| [1 0 0 -1 0 1 1 0] [r3] |
| [0 0 0 0 1 -1 1 -1] * [i0] |
| [0 -1 1 0 1 0 0 1] [i1] |
| [0 0 0 0 1 1 -1 -1] [i2] |
| [0 1 -1 0 1 0 0 1] [i3] |
| */ |
| |
| v4sf sr0 = VADD(r0,r3), dr0 = VSUB(r0,r3); |
| v4sf sr1 = VADD(r1,r2), dr1 = VSUB(r1,r2); |
| v4sf si0 = VADD(i0,i3), di0 = VSUB(i0,i3); |
| v4sf si1 = VADD(i1,i2), di1 = VSUB(i1,i2); |
| |
| r0 = VADD(sr0, sr1); |
| r2 = VSUB(sr0, sr1); |
| r1 = VSUB(dr0, si1); |
| r3 = VADD(dr0, si1); |
| i0 = VSUB(di0, di1); |
| i2 = VADD(di0, di1); |
| i1 = VSUB(si0, dr1); |
| i3 = VADD(si0, dr1); |
| |
| VCPLXMULCONJ(r1,i1,e[0],e[1]); |
| VCPLXMULCONJ(r2,i2,e[2],e[3]); |
| VCPLXMULCONJ(r3,i3,e[4],e[5]); |
| |
| VTRANSPOSE4(r0,r1,r2,r3); |
| VTRANSPOSE4(i0,i1,i2,i3); |
| |
| if (!first) { |
| *out++ = r0; |
| *out++ = i0; |
| } |
| *out++ = r1; |
| *out++ = i1; |
| *out++ = r2; |
| *out++ = i2; |
| *out++ = r3; |
| *out++ = i3; |
| } |
| |
| static NEVER_INLINE(void) pffft_real_preprocess(int Ncvec, const v4sf *in, v4sf *out, const v4sf *e) { |
| int k, dk = Ncvec/SIMD_SZ; // number of 4x4 matrix blocks |
| /* fftpack order is f0r f1r f1i f2r f2i ... f(n-1)r f(n-1)i f(n)r */ |
| |
| v4sf_union Xr, Xi, *uout = (v4sf_union*)out; |
| float cr0, ci0, cr1, ci1, cr2, ci2, cr3, ci3; |
| static const float s = (float)M_SQRT2; |
| assert(in != out); |
| for (k=0; k < 4; ++k) { |
| Xr.f[k] = ((float*)in)[8*k]; |
| Xi.f[k] = ((float*)in)[8*k+4]; |
| } |
| |
| pffft_real_preprocess_4x4(in, e, out+1, 1); // will write only 6 values |
| |
| /* |
| [Xr0 Xr1 Xr2 Xr3 Xi0 Xi1 Xi2 Xi3] |
| |
| [cr0] [1 0 2 0 1 0 0 0] |
| [cr1] [1 0 0 0 -1 0 -2 0] |
| [cr2] [1 0 -2 0 1 0 0 0] |
| [cr3] [1 0 0 0 -1 0 2 0] |
| [ci0] [0 2 0 2 0 0 0 0] |
| [ci1] [0 s 0 -s 0 -s 0 -s] |
| [ci2] [0 0 0 0 0 -2 0 2] |
| [ci3] [0 -s 0 s 0 -s 0 -s] |
| */ |
| for (k=1; k < dk; ++k) { |
| pffft_real_preprocess_4x4(in+8*k, e + k*6, out-1+k*8, 0); |
| } |
| |
| cr0=(Xr.f[0]+Xi.f[0]) + 2*Xr.f[2]; uout[0].f[0] = cr0; |
| cr1=(Xr.f[0]-Xi.f[0]) - 2*Xi.f[2]; uout[0].f[1] = cr1; |
| cr2=(Xr.f[0]+Xi.f[0]) - 2*Xr.f[2]; uout[0].f[2] = cr2; |
| cr3=(Xr.f[0]-Xi.f[0]) + 2*Xi.f[2]; uout[0].f[3] = cr3; |
| ci0= 2*(Xr.f[1]+Xr.f[3]); uout[2*Ncvec-1].f[0] = ci0; |
| ci1= s*(Xr.f[1]-Xr.f[3]) - s*(Xi.f[1]+Xi.f[3]); uout[2*Ncvec-1].f[1] = ci1; |
| ci2= 2*(Xi.f[3]-Xi.f[1]); uout[2*Ncvec-1].f[2] = ci2; |
| ci3=-s*(Xr.f[1]-Xr.f[3]) - s*(Xi.f[1]+Xi.f[3]); uout[2*Ncvec-1].f[3] = ci3; |
| } |
| |
| |
| void pffft_transform_internal(PFFFT_Setup *setup, const float *finput, float *foutput, v4sf *scratch, |
| pffft_direction_t direction, int ordered) { |
| int k, Ncvec = setup->Ncvec; |
| int nf_odd = (setup->ifac[1] & 1); |
| |
| // temporary buffer is allocated on the stack if the scratch pointer is NULL |
| int stack_allocate = (scratch == 0 ? Ncvec*2 : 1); |
| VLA_ARRAY_ON_STACK(v4sf, scratch_on_stack, stack_allocate); |
| |
| const v4sf *vinput = (const v4sf*)finput; |
| v4sf *voutput = (v4sf*)foutput; |
| v4sf *buff[2] = { voutput, scratch ? scratch : scratch_on_stack }; |
| int ib = (nf_odd ^ ordered ? 1 : 0); |
| |
| assert(VALIGNED(finput) && VALIGNED(foutput)); |
| |
| //assert(finput != foutput); |
| if (direction == PFFFT_FORWARD) { |
| ib = !ib; |
| if (setup->transform == PFFFT_REAL) { |
| ib = (rfftf1_ps(Ncvec*2, vinput, buff[ib], buff[!ib], |
| setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1); |
| pffft_real_finalize(Ncvec, buff[ib], buff[!ib], (v4sf*)setup->e); |
| } else { |
| v4sf *tmp = buff[ib]; |
| for (k=0; k < Ncvec; ++k) { |
| UNINTERLEAVE2(vinput[k*2], vinput[k*2+1], tmp[k*2], tmp[k*2+1]); |
| } |
| ib = (cfftf1_ps(Ncvec, buff[ib], buff[!ib], buff[ib], |
| setup->twiddle, &setup->ifac[0], -1) == buff[0] ? 0 : 1); |
| pffft_cplx_finalize(Ncvec, buff[ib], buff[!ib], (v4sf*)setup->e); |
| } |
| if (ordered) { |
| pffft_zreorder(setup, (float*)buff[!ib], (float*)buff[ib], PFFFT_FORWARD); |
| } else ib = !ib; |
| } else { |
| if (vinput == buff[ib]) { |
| ib = !ib; // may happen when finput == foutput |
| } |
| if (ordered) { |
| pffft_zreorder(setup, (float*)vinput, (float*)buff[ib], PFFFT_BACKWARD); |
| vinput = buff[ib]; ib = !ib; |
| } |
| if (setup->transform == PFFFT_REAL) { |
| pffft_real_preprocess(Ncvec, vinput, buff[ib], (v4sf*)setup->e); |
| ib = (rfftb1_ps(Ncvec*2, buff[ib], buff[0], buff[1], |
| setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1); |
| } else { |
| pffft_cplx_preprocess(Ncvec, vinput, buff[ib], (v4sf*)setup->e); |
| ib = (cfftf1_ps(Ncvec, buff[ib], buff[0], buff[1], |
| setup->twiddle, &setup->ifac[0], +1) == buff[0] ? 0 : 1); |
| for (k=0; k < Ncvec; ++k) { |
| INTERLEAVE2(buff[ib][k*2], buff[ib][k*2+1], buff[ib][k*2], buff[ib][k*2+1]); |
| } |
| } |
| } |
| |
| if (buff[ib] != voutput) { |
| /* extra copy required -- this situation should only happen when finput == foutput */ |
| assert(finput==foutput); |
| for (k=0; k < Ncvec; ++k) { |
| v4sf a = buff[ib][2*k], b = buff[ib][2*k+1]; |
| voutput[2*k] = a; voutput[2*k+1] = b; |
| } |
| ib = !ib; |
| } |
| assert(buff[ib] == voutput); |
| } |
| |
| void pffft_zconvolve_accumulate(PFFFT_Setup *s, const float *a, const float *b, float *ab, float scaling) { |
| int Ncvec = s->Ncvec; |
| const v4sf * RESTRICT va = (const v4sf*)a; |
| const v4sf * RESTRICT vb = (const v4sf*)b; |
| v4sf * RESTRICT vab = (v4sf*)ab; |
| |
| #ifdef __arm__ |
| __builtin_prefetch(va); |
| __builtin_prefetch(vb); |
| __builtin_prefetch(vab); |
| __builtin_prefetch(va+2); |
| __builtin_prefetch(vb+2); |
| __builtin_prefetch(vab+2); |
| __builtin_prefetch(va+4); |
| __builtin_prefetch(vb+4); |
| __builtin_prefetch(vab+4); |
| __builtin_prefetch(va+6); |
| __builtin_prefetch(vb+6); |
| __builtin_prefetch(vab+6); |
| # ifndef __clang__ |
| # define ZCONVOLVE_USING_INLINE_NEON_ASM |
| # endif |
| #endif |
| |
| float ar, ai, br, bi, abr, abi; |
| #ifndef ZCONVOLVE_USING_INLINE_ASM |
| v4sf vscal = LD_PS1(scaling); |
| int i; |
| #endif |
| |
| assert(VALIGNED(a) && VALIGNED(b) && VALIGNED(ab)); |
| ar = ((v4sf_union*)va)[0].f[0]; |
| ai = ((v4sf_union*)va)[1].f[0]; |
| br = ((v4sf_union*)vb)[0].f[0]; |
| bi = ((v4sf_union*)vb)[1].f[0]; |
| abr = ((v4sf_union*)vab)[0].f[0]; |
| abi = ((v4sf_union*)vab)[1].f[0]; |
| |
| #ifdef ZCONVOLVE_USING_INLINE_ASM // inline asm version, unfortunately miscompiled by clang 3.2, at least on ubuntu.. so this will be restricted to gcc |
| const float *a_ = a, *b_ = b; float *ab_ = ab; |
| int N = Ncvec; |
| asm volatile("mov r8, %2 \n" |
| "vdup.f32 q15, %4 \n" |
| "1: \n" |
| "pld [%0,#64] \n" |
| "pld [%1,#64] \n" |
| "pld [%2,#64] \n" |
| "pld [%0,#96] \n" |
| "pld [%1,#96] \n" |
| "pld [%2,#96] \n" |
| "vld1.f32 {q0,q1}, [%0,:128]! \n" |
| "vld1.f32 {q4,q5}, [%1,:128]! \n" |
| "vld1.f32 {q2,q3}, [%0,:128]! \n" |
| "vld1.f32 {q6,q7}, [%1,:128]! \n" |
| "vld1.f32 {q8,q9}, [r8,:128]! \n" |
| |
| "vmul.f32 q10, q0, q4 \n" |
| "vmul.f32 q11, q0, q5 \n" |
| "vmul.f32 q12, q2, q6 \n" |
| "vmul.f32 q13, q2, q7 \n" |
| "vmls.f32 q10, q1, q5 \n" |
| "vmla.f32 q11, q1, q4 \n" |
| "vld1.f32 {q0,q1}, [r8,:128]! \n" |
| "vmls.f32 q12, q3, q7 \n" |
| "vmla.f32 q13, q3, q6 \n" |
| "vmla.f32 q8, q10, q15 \n" |
| "vmla.f32 q9, q11, q15 \n" |
| "vmla.f32 q0, q12, q15 \n" |
| "vmla.f32 q1, q13, q15 \n" |
| "vst1.f32 {q8,q9},[%2,:128]! \n" |
| "vst1.f32 {q0,q1},[%2,:128]! \n" |
| "subs %3, #2 \n" |
| "bne 1b \n" |
| : "+r"(a_), "+r"(b_), "+r"(ab_), "+r"(N) : "r"(scaling) : "r8", "q0","q1","q2","q3","q4","q5","q6","q7","q8","q9", "q10","q11","q12","q13","q15","memory"); |
| #else // default routine, works fine for non-arm cpus with current compilers |
| for (i=0; i < Ncvec; i += 2) { |
| v4sf ar, ai, br, bi; |
| ar = va[2*i+0]; ai = va[2*i+1]; |
| br = vb[2*i+0]; bi = vb[2*i+1]; |
| VCPLXMUL(ar, ai, br, bi); |
| vab[2*i+0] = VMADD(ar, vscal, vab[2*i+0]); |
| vab[2*i+1] = VMADD(ai, vscal, vab[2*i+1]); |
| ar = va[2*i+2]; ai = va[2*i+3]; |
| br = vb[2*i+2]; bi = vb[2*i+3]; |
| VCPLXMUL(ar, ai, br, bi); |
| vab[2*i+2] = VMADD(ar, vscal, vab[2*i+2]); |
| vab[2*i+3] = VMADD(ai, vscal, vab[2*i+3]); |
| } |
| #endif |
| if (s->transform == PFFFT_REAL) { |
| ((v4sf_union*)vab)[0].f[0] = abr + ar*br*scaling; |
| ((v4sf_union*)vab)[1].f[0] = abi + ai*bi*scaling; |
| } |
| } |
| |
| |
| void pffft_zconvolve_no_accu(PFFFT_Setup *s, const float *a, const float *b, float *ab, float scaling) { |
| v4sf vscal = LD_PS1(scaling); |
| const v4sf * RESTRICT va = (const v4sf*)a; |
| const v4sf * RESTRICT vb = (const v4sf*)b; |
| v4sf * RESTRICT vab = (v4sf*)ab; |
| float sar, sai, sbr, sbi; |
| const int NcvecMulTwo = 2*s->Ncvec; /* int Ncvec = s->Ncvec; */ |
| int k; /* was i -- but always used "2*i" - except at for() */ |
| |
| #ifdef __arm__ |
| __builtin_prefetch(va); |
| __builtin_prefetch(vb); |
| __builtin_prefetch(vab); |
| __builtin_prefetch(va+2); |
| __builtin_prefetch(vb+2); |
| __builtin_prefetch(vab+2); |
| __builtin_prefetch(va+4); |
| __builtin_prefetch(vb+4); |
| __builtin_prefetch(vab+4); |
| __builtin_prefetch(va+6); |
| __builtin_prefetch(vb+6); |
| __builtin_prefetch(vab+6); |
| # ifndef __clang__ |
| # define ZCONVOLVE_USING_INLINE_NEON_ASM |
| # endif |
| #endif |
| |
| assert(VALIGNED(a) && VALIGNED(b) && VALIGNED(ab)); |
| sar = ((v4sf_union*)va)[0].f[0]; |
| sai = ((v4sf_union*)va)[1].f[0]; |
| sbr = ((v4sf_union*)vb)[0].f[0]; |
| sbi = ((v4sf_union*)vb)[1].f[0]; |
| |
| /* default routine, works fine for non-arm cpus with current compilers */ |
| for (k=0; k < NcvecMulTwo; k += 4) { |
| v4sf var, vai, vbr, vbi; |
| var = va[k+0]; vai = va[k+1]; |
| vbr = vb[k+0]; vbi = vb[k+1]; |
| VCPLXMUL(var, vai, vbr, vbi); |
| vab[k+0] = VMUL(var, vscal); |
| vab[k+1] = VMUL(vai, vscal); |
| var = va[k+2]; vai = va[k+3]; |
| vbr = vb[k+2]; vbi = vb[k+3]; |
| VCPLXMUL(var, vai, vbr, vbi); |
| vab[k+2] = VMUL(var, vscal); |
| vab[k+3] = VMUL(vai, vscal); |
| } |
| |
| if (s->transform == PFFFT_REAL) { |
| ((v4sf_union*)vab)[0].f[0] = sar*sbr*scaling; |
| ((v4sf_union*)vab)[1].f[0] = sai*sbi*scaling; |
| } |
| } |
| |
| |
| #else // defined(PFFFT_SIMD_DISABLE) |
| |
| // standard routine using scalar floats, without SIMD stuff. |
| |
| #define pffft_zreorder_nosimd pffft_zreorder |
| void pffft_zreorder_nosimd(PFFFT_Setup *setup, const float *in, float *out, pffft_direction_t direction) { |
| int k, N = setup->N; |
| if (setup->transform == PFFFT_COMPLEX) { |
| for (k=0; k < 2*N; ++k) out[k] = in[k]; |
| return; |
| } |
| else if (direction == PFFFT_FORWARD) { |
| float x_N = in[N-1]; |
| for (k=N-1; k > 1; --k) out[k] = in[k-1]; |
| out[0] = in[0]; |
| out[1] = x_N; |
| } else { |
| float x_N = in[1]; |
| for (k=1; k < N-1; ++k) out[k] = in[k+1]; |
| out[0] = in[0]; |
| out[N-1] = x_N; |
| } |
| } |
| |
| #define pffft_transform_internal_nosimd pffft_transform_internal |
| void pffft_transform_internal_nosimd(PFFFT_Setup *setup, const float *input, float *output, float *scratch, |
| pffft_direction_t direction, int ordered) { |
| int Ncvec = setup->Ncvec; |
| int nf_odd = (setup->ifac[1] & 1); |
| |
| // temporary buffer is allocated on the stack if the scratch pointer is NULL |
| int stack_allocate = (scratch == 0 ? Ncvec*2 : 1); |
| VLA_ARRAY_ON_STACK(v4sf, scratch_on_stack, stack_allocate); |
| float *buff[2]; |
| int ib; |
| if (scratch == 0) scratch = scratch_on_stack; |
| buff[0] = output; buff[1] = scratch; |
| |
| if (setup->transform == PFFFT_COMPLEX) ordered = 0; // it is always ordered. |
| ib = (nf_odd ^ ordered ? 1 : 0); |
| |
| if (direction == PFFFT_FORWARD) { |
| if (setup->transform == PFFFT_REAL) { |
| ib = (rfftf1_ps(Ncvec*2, input, buff[ib], buff[!ib], |
| setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1); |
| } else { |
| ib = (cfftf1_ps(Ncvec, input, buff[ib], buff[!ib], |
| setup->twiddle, &setup->ifac[0], -1) == buff[0] ? 0 : 1); |
| } |
| if (ordered) { |
| pffft_zreorder(setup, buff[ib], buff[!ib], PFFFT_FORWARD); ib = !ib; |
| } |
| } else { |
| if (input == buff[ib]) { |
| ib = !ib; // may happen when finput == foutput |
| } |
| if (ordered) { |
| pffft_zreorder(setup, input, buff[!ib], PFFFT_BACKWARD); |
| input = buff[!ib]; |
| } |
| if (setup->transform == PFFFT_REAL) { |
| ib = (rfftb1_ps(Ncvec*2, input, buff[ib], buff[!ib], |
| setup->twiddle, &setup->ifac[0]) == buff[0] ? 0 : 1); |
| } else { |
| ib = (cfftf1_ps(Ncvec, input, buff[ib], buff[!ib], |
| setup->twiddle, &setup->ifac[0], +1) == buff[0] ? 0 : 1); |
| } |
| } |
| if (buff[ib] != output) { |
| int k; |
| // extra copy required -- this situation should happens only when finput == foutput |
| assert(input==output); |
| for (k=0; k < Ncvec; ++k) { |
| float a = buff[ib][2*k], b = buff[ib][2*k+1]; |
| output[2*k] = a; output[2*k+1] = b; |
| } |
| ib = !ib; |
| } |
| assert(buff[ib] == output); |
| } |
| |
| #define pffft_zconvolve_accumulate_nosimd pffft_zconvolve_accumulate |
| void pffft_zconvolve_accumulate_nosimd(PFFFT_Setup *s, const float *a, const float *b, |
| float *ab, float scaling) { |
| int NcvecMulTwo = 2*s->Ncvec; /* int Ncvec = s->Ncvec; */ |
| int k; /* was i -- but always used "2*i" - except at for() */ |
| |
| if (s->transform == PFFFT_REAL) { |
| // take care of the fftpack ordering |
| ab[0] += a[0]*b[0]*scaling; |
| ab[NcvecMulTwo-1] += a[NcvecMulTwo-1]*b[NcvecMulTwo-1]*scaling; |
| ++ab; ++a; ++b; NcvecMulTwo -= 2; |
| } |
| for (k=0; k < NcvecMulTwo; k += 2) { |
| float ar, ai, br, bi; |
| ar = a[k+0]; ai = a[k+1]; |
| br = b[k+0]; bi = b[k+1]; |
| VCPLXMUL(ar, ai, br, bi); |
| ab[k+0] += ar*scaling; |
| ab[k+1] += ai*scaling; |
| } |
| } |
| |
| |
| #define pffft_zconvolve_no_accu_nosimd pffft_zconvolve_no_accu |
| void pffft_zconvolve_no_accu_nosimd(PFFFT_Setup *s, const float *a, const float *b, |
| float *ab, float scaling) { |
| int NcvecMulTwo = 2*s->Ncvec; /* int Ncvec = s->Ncvec; */ |
| int k; /* was i -- but always used "2*i" - except at for() */ |
| |
| if (s->transform == PFFFT_REAL) { |
| // take care of the fftpack ordering |
| ab[0] += a[0]*b[0]*scaling; |
| ab[NcvecMulTwo-1] += a[NcvecMulTwo-1]*b[NcvecMulTwo-1]*scaling; |
| ++ab; ++a; ++b; NcvecMulTwo -= 2; |
| } |
| for (k=0; k < NcvecMulTwo; k += 2) { |
| float ar, ai, br, bi; |
| ar = a[k+0]; ai = a[k+1]; |
| br = b[k+0]; bi = b[k+1]; |
| VCPLXMUL(ar, ai, br, bi); |
| ab[k+0] = ar*scaling; |
| ab[k+1] = ai*scaling; |
| } |
| } |
| |
| |
| #endif // defined(PFFFT_SIMD_DISABLE) |
| |
| void pffft_transform(PFFFT_Setup *setup, const float *input, float *output, float *work, pffft_direction_t direction) { |
| pffft_transform_internal(setup, input, output, (v4sf*)work, direction, 0); |
| } |
| |
| void pffft_transform_ordered(PFFFT_Setup *setup, const float *input, float *output, float *work, pffft_direction_t direction) { |
| pffft_transform_internal(setup, input, output, (v4sf*)work, direction, 1); |
| } |
| |
| |
| int pffft_next_power_of_two(int N) { |
| /* https://graphics.stanford.edu/~seander/bithacks.html#RoundUpPowerOf2 */ |
| /* compute the next highest power of 2 of 32-bit v */ |
| unsigned v = N; |
| v--; |
| v |= v >> 1; |
| v |= v >> 2; |
| v |= v >> 4; |
| v |= v >> 8; |
| v |= v >> 16; |
| v++; |
| return v; |
| } |
| |
| int pffft_is_power_of_two(int N) { |
| /* https://graphics.stanford.edu/~seander/bithacks.html#DetermineIfPowerOf2 */ |
| int f = N && !(N & (N - 1)); |
| return f; |
| } |
| |
| int pffft_min_fft_size(pffft_transform_t transform) { |
| /* unfortunately, the fft size must be a multiple of 16 for complex FFTs |
| and 32 for real FFTs -- a lot of stuff would need to be rewritten to |
| handle other cases (or maybe just switch to a scalar fft, I don't know..) */ |
| if (transform == PFFFT_REAL) |
| return ( 2 * SIMD_SZ * SIMD_SZ ); |
| else if (transform == PFFFT_COMPLEX) |
| return ( SIMD_SZ * SIMD_SZ ); |
| else |
| return 1; |
| } |
| |