Add the Ooura FFT to RealFourier.

We are using the Ooura FFT in a few places:
- AGC
- Transient suppression
- Noise suppression

The optimized OpenMAX DL FFT is considerably faster, but currently does
not compile everywhere, notably on iOS. This change will allow us to use
Openmax when possible and otherwise fall back to Ooura.

(Unfortunately, noise suppression won't be able to take advantage of it
since it's not C++. Upgrade time?)

R=aluebs@webrtc.org, mgraczyk@chromium.org

Review URL: https://webrtc-codereview.appspot.com/45789004

Cr-Commit-Position: refs/heads/master@{#8798}
git-svn-id: http://webrtc.googlecode.com/svn/trunk@8798 4adac7df-926f-26a2-2b94-8c16560cd09d

7 files changed, 6 insertions, 1353 deletions

diff --git a/webrtc/modules/audio_processing/BUILD.gn b/webrtc/modules/audio_processing/BUILD.gn
index b34c41cc32..e4ee78266c 100644
--- a/webrtc/modules/audio_processing/BUILD.gn
+++ b/webrtc/modules/audio_processing/BUILD.gn
@@ -116,8 +116,6 @@ source_set("audio_processing") {
     "utility/delay_estimator_internal.h",
     "utility/delay_estimator_wrapper.c",
     "utility/delay_estimator_wrapper.h",
-    "utility/fft4g.c",
-    "utility/fft4g.h",
     "voice_detection_impl.cc",
     "voice_detection_impl.h",
   ]
diff --git a/webrtc/modules/audio_processing/agc/agc_audio_proc.cc b/webrtc/modules/audio_processing/agc/agc_audio_proc.cc
index ccc77cc9ea..dc4a5a711c 100644
--- a/webrtc/modules/audio_processing/agc/agc_audio_proc.cc
+++ b/webrtc/modules/audio_processing/agc/agc_audio_proc.cc
@@ -13,6 +13,7 @@
 #include <math.h>
 #include <stdio.h>
 
+#include "webrtc/common_audio/fft4g.h"
 #include "webrtc/modules/audio_processing/agc/agc_audio_proc_internal.h"
 #include "webrtc/modules/audio_processing/agc/pitch_internal.h"
 #include "webrtc/modules/audio_processing/agc/pole_zero_filter.h"
@@ -21,7 +22,6 @@ extern "C" {
 #include "webrtc/modules/audio_coding/codecs/isac/main/source/lpc_analysis.h"
 #include "webrtc/modules/audio_coding/codecs/isac/main/source/pitch_estimator.h"
 #include "webrtc/modules/audio_coding/codecs/isac/main/source/structs.h"
-#include "webrtc/modules/audio_processing/utility/fft4g.h"
 }
 #include "webrtc/modules/interface/module_common_types.h"
 
diff --git a/webrtc/modules/audio_processing/audio_processing.gypi b/webrtc/modules/audio_processing/audio_processing.gypi
index 3175cc66e4..9d1ac8ccbd 100644
--- a/webrtc/modules/audio_processing/audio_processing.gypi
+++ b/webrtc/modules/audio_processing/audio_processing.gypi
@@ -126,8 +126,6 @@
         'utility/delay_estimator_internal.h',
         'utility/delay_estimator_wrapper.c',
         'utility/delay_estimator_wrapper.h',
-        'utility/fft4g.c',
-        'utility/fft4g.h',
         'voice_detection_impl.cc',
         'voice_detection_impl.h',
       ],
diff --git a/webrtc/modules/audio_processing/ns/ns_core.c b/webrtc/modules/audio_processing/ns/ns_core.c
index 56f2f46aa6..9e230dd140 100644
--- a/webrtc/modules/audio_processing/ns/ns_core.c
+++ b/webrtc/modules/audio_processing/ns/ns_core.c
@@ -13,11 +13,11 @@
 #include <string.h>
 #include <stdlib.h>
 
+#include "webrtc/common_audio/fft4g.h"
 #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
 #include "webrtc/modules/audio_processing/ns/include/noise_suppression.h"
 #include "webrtc/modules/audio_processing/ns/ns_core.h"
 #include "webrtc/modules/audio_processing/ns/windows_private.h"
-#include "webrtc/modules/audio_processing/utility/fft4g.h"
 
 // Set Feature Extraction Parameters.
 static void set_feature_extraction_parameters(NoiseSuppressionC* self) {
diff --git a/webrtc/modules/audio_processing/transient/transient_suppressor.cc b/webrtc/modules/audio_processing/transient/transient_suppressor.cc
index 44ea6bf233..2f79a20ac7 100644
--- a/webrtc/modules/audio_processing/transient/transient_suppressor.cc
+++ b/webrtc/modules/audio_processing/transient/transient_suppressor.cc
@@ -18,14 +18,12 @@
 #include <set>
 
 #include "webrtc/base/scoped_ptr.h"
+#include "webrtc/common_audio/fft4g.h"
 #include "webrtc/common_audio/include/audio_util.h"
 #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
 #include "webrtc/modules/audio_processing/transient/common.h"
 #include "webrtc/modules/audio_processing/transient/transient_detector.h"
 #include "webrtc/modules/audio_processing/ns/windows_private.h"
-extern "C" {
-#include "webrtc/modules/audio_processing/utility/fft4g.h"
-}
 #include "webrtc/system_wrappers/interface/logging.h"
 #include "webrtc/typedefs.h"
 
@@ -39,10 +37,12 @@ static const size_t kMinVoiceBin = 3;
 static const size_t kMaxVoiceBin = 60;
 
 namespace {
+
 float ComplexMagnitude(float a, float b) {
   return std::abs(a) + std::abs(b);
 }
-}
+
+}  // namespace
 
 TransientSuppressor::TransientSuppressor()
     : data_length_(0),
diff --git a/webrtc/modules/audio_processing/utility/fft4g.c b/webrtc/modules/audio_processing/utility/fft4g.c
deleted file mode 100644
index cbc4dc31eb..0000000000
--- a/webrtc/modules/audio_processing/utility/fft4g.c
+++ /dev/null
@@ -1,1326 +0,0 @@
-/*
- * http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html
- * Copyright Takuya OOURA, 1996-2001
- *
- * You may use, copy, modify and distribute this code for any purpose (include
- * commercial use) and without fee. Please refer to this package when you modify
- * this code.
- *
- * Changes:
- * Trivial type modifications by the WebRTC authors.
- */
-
-/*
-Fast Fourier/Cosine/Sine Transform
-    dimension   :one
-    data length :power of 2
-    decimation  :frequency
-    radix       :4, 2
-    data        :inplace
-    table       :use
-functions
-    cdft: Complex Discrete Fourier Transform
-    rdft: Real Discrete Fourier Transform
-    ddct: Discrete Cosine Transform
-    ddst: Discrete Sine Transform
-    dfct: Cosine Transform of RDFT (Real Symmetric DFT)
-    dfst: Sine Transform of RDFT (Real Anti-symmetric DFT)
-function prototypes
-    void cdft(int, int, float *, int *, float *);
-    void rdft(int, int, float *, int *, float *);
-    void ddct(int, int, float *, int *, float *);
-    void ddst(int, int, float *, int *, float *);
-    void dfct(int, float *, float *, int *, float *);
-    void dfst(int, float *, float *, int *, float *);
-
-
--------- Complex DFT (Discrete Fourier Transform) --------
-    [definition]
-        <case1>
-            X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
-        <case2>
-            X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
-        (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            cdft(2*n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            cdft(2*n, -1, a, ip, w);
-    [parameters]
-        2*n            :data length (int)
-                        n >= 1, n = power of 2
-        a[0...2*n-1]   :input/output data (float *)
-                        input data
-                            a[2*j] = Re(x[j]),
-                            a[2*j+1] = Im(x[j]), 0<=j<n
-                        output data
-                            a[2*k] = Re(X[k]),
-                            a[2*k+1] = Im(X[k]), 0<=k<n
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n)
-                        strictly,
-                        length of ip >=
-                            2+(1<<(int)(log(n+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n/2-1]   :cos/sin table (float *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of
-            cdft(2*n, -1, a, ip, w);
-        is
-            cdft(2*n, 1, a, ip, w);
-            for (j = 0; j <= 2 * n - 1; j++) {
-                a[j] *= 1.0 / n;
-            }
-        .
-
-
--------- Real DFT / Inverse of Real DFT --------
-    [definition]
-        <case1> RDFT
-            R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
-            I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
-        <case2> IRDFT (excluding scale)
-            a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
-                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
-                   sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            rdft(n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            rdft(n, -1, a, ip, w);
-    [parameters]
-        n              :data length (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (float *)
-                        <case1>
-                            output data
-                                a[2*k] = R[k], 0<=k<n/2
-                                a[2*k+1] = I[k], 0<k<n/2
-                                a[1] = R[n/2]
-                        <case2>
-                            input data
-                                a[2*j] = R[j], 0<=j<n/2
-                                a[2*j+1] = I[j], 0<j<n/2
-                                a[1] = R[n/2]
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/2)
-                        strictly,
-                        length of ip >=
-                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n/2-1]   :cos/sin table (float *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of
-            rdft(n, 1, a, ip, w);
-        is
-            rdft(n, -1, a, ip, w);
-            for (j = 0; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- DCT (Discrete Cosine Transform) / Inverse of DCT --------
-    [definition]
-        <case1> IDCT (excluding scale)
-            C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
-        <case2> DCT
-            C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            ddct(n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            ddct(n, -1, a, ip, w);
-    [parameters]
-        n              :data length (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (float *)
-                        output data
-                            a[k] = C[k], 0<=k<n
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/2)
-                        strictly,
-                        length of ip >=
-                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/4-1] :cos/sin table (float *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of
-            ddct(n, -1, a, ip, w);
-        is
-            a[0] *= 0.5;
-            ddct(n, 1, a, ip, w);
-            for (j = 0; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- DST (Discrete Sine Transform) / Inverse of DST --------
-    [definition]
-        <case1> IDST (excluding scale)
-            S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
-        <case2> DST
-            S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
-    [usage]
-        <case1>
-            ip[0] = 0; // first time only
-            ddst(n, 1, a, ip, w);
-        <case2>
-            ip[0] = 0; // first time only
-            ddst(n, -1, a, ip, w);
-    [parameters]
-        n              :data length (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (float *)
-                        <case1>
-                            input data
-                                a[j] = A[j], 0<j<n
-                                a[0] = A[n]
-                            output data
-                                a[k] = S[k], 0<=k<n
-                        <case2>
-                            output data
-                                a[k] = S[k], 0<k<n
-                                a[0] = S[n]
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/2)
-                        strictly,
-                        length of ip >=
-                            2+(1<<(int)(log(n/2+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/4-1] :cos/sin table (float *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of
-            ddst(n, -1, a, ip, w);
-        is
-            a[0] *= 0.5;
-            ddst(n, 1, a, ip, w);
-            for (j = 0; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- Cosine Transform of RDFT (Real Symmetric DFT) --------
-    [definition]
-        C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n
-    [usage]
-        ip[0] = 0; // first time only
-        dfct(n, a, t, ip, w);
-    [parameters]
-        n              :data length - 1 (int)
-                        n >= 2, n = power of 2
-        a[0...n]       :input/output data (float *)
-                        output data
-                            a[k] = C[k], 0<=k<=n
-        t[0...n/2]     :work area (float *)
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/4)
-                        strictly,
-                        length of ip >=
-                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/8-1] :cos/sin table (float *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of
-            a[0] *= 0.5;
-            a[n] *= 0.5;
-            dfct(n, a, t, ip, w);
-        is
-            a[0] *= 0.5;
-            a[n] *= 0.5;
-            dfct(n, a, t, ip, w);
-            for (j = 0; j <= n; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
--------- Sine Transform of RDFT (Real Anti-symmetric DFT) --------
-    [definition]
-        S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n
-    [usage]
-        ip[0] = 0; // first time only
-        dfst(n, a, t, ip, w);
-    [parameters]
-        n              :data length + 1 (int)
-                        n >= 2, n = power of 2
-        a[0...n-1]     :input/output data (float *)
-                        output data
-                            a[k] = S[k], 0<k<n
-                        (a[0] is used for work area)
-        t[0...n/2-1]   :work area (float *)
-        ip[0...*]      :work area for bit reversal (int *)
-                        length of ip >= 2+sqrt(n/4)
-                        strictly,
-                        length of ip >=
-                            2+(1<<(int)(log(n/4+0.5)/log(2))/2).
-                        ip[0],ip[1] are pointers of the cos/sin table.
-        w[0...n*5/8-1] :cos/sin table (float *)
-                        w[],ip[] are initialized if ip[0] == 0.
-    [remark]
-        Inverse of
-            dfst(n, a, t, ip, w);
-        is
-            dfst(n, a, t, ip, w);
-            for (j = 1; j <= n - 1; j++) {
-                a[j] *= 2.0 / n;
-            }
-        .
-
-
-Appendix :
-    The cos/sin table is recalculated when the larger table required.
-    w[] and ip[] are compatible with all routines.
-*/
-
-static void makewt(int nw, int *ip, float *w);
-static void makect(int nc, int *ip, float *c);
-static void bitrv2(int n, int *ip, float *a);
-static void bitrv2conj(int n, int *ip, float *a);
-static void cftfsub(int n, float *a, float *w);
-static void cftbsub(int n, float *a, float *w);
-static void cft1st(int n, float *a, float *w);
-static void cftmdl(int n, int l, float *a, float *w);
-static void rftfsub(int n, float *a, int nc, float *c);
-static void rftbsub(int n, float *a, int nc, float *c);
-#if 0  // Not used.
-static void dctsub(int n, float *a, int nc, float *c)
-static void dstsub(int n, float *a, int nc, float *c)
-#endif
-
-
-void WebRtc_cdft(int n, int isgn, float *a, int *ip, float *w)
-{
-    if (n > (ip[0] << 2)) {
-        makewt(n >> 2, ip, w);
-    }
-    if (n > 4) {
-        if (isgn >= 0) {
-            bitrv2(n, ip + 2, a);
-            cftfsub(n, a, w);
-        } else {
-            bitrv2conj(n, ip + 2, a);
-            cftbsub(n, a, w);
-        }
-    } else if (n == 4) {
-        cftfsub(n, a, w);
-    }
-}
-
-
-void WebRtc_rdft(int n, int isgn, float *a, int *ip, float *w)
-{
-    int nw, nc;
-    float xi;
-
-    nw = ip[0];
-    if (n > (nw << 2)) {
-        nw = n >> 2;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > (nc << 2)) {
-        nc = n >> 2;
-        makect(nc, ip, w + nw);
-    }
-    if (isgn >= 0) {
-        if (n > 4) {
-            bitrv2(n, ip + 2, a);
-            cftfsub(n, a, w);
-            rftfsub(n, a, nc, w + nw);
-        } else if (n == 4) {
-            cftfsub(n, a, w);
-        }
-        xi = a[0] - a[1];
-        a[0] += a[1];
-        a[1] = xi;
-    } else {
-        a[1] = 0.5f * (a[0] - a[1]);
-        a[0] -= a[1];
-        if (n > 4) {
-            rftbsub(n, a, nc, w + nw);
-            bitrv2(n, ip + 2, a);
-            cftbsub(n, a, w);
-        } else if (n == 4) {
-            cftfsub(n, a, w);
-        }
-    }
-}
-
-#if 0  // Not used.
-static void ddct(int n, int isgn, float *a, int *ip, float *w)
-{
-    int j, nw, nc;
-    float xr;
-
-    nw = ip[0];
-    if (n > (nw << 2)) {
-        nw = n >> 2;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > nc) {
-        nc = n;
-        makect(nc, ip, w + nw);
-    }
-    if (isgn < 0) {
-        xr = a[n - 1];
-        for (j = n - 2; j >= 2; j -= 2) {
-            a[j + 1] = a[j] - a[j - 1];
-            a[j] += a[j - 1];
-        }
-        a[1] = a[0] - xr;
-        a[0] += xr;
-        if (n > 4) {
-            rftbsub(n, a, nc, w + nw);
-            bitrv2(n, ip + 2, a);
-            cftbsub(n, a, w);
-        } else if (n == 4) {
-            cftfsub(n, a, w);
-        }
-    }
-    dctsub(n, a, nc, w + nw);
-    if (isgn >= 0) {
-        if (n > 4) {
-            bitrv2(n, ip + 2, a);
-            cftfsub(n, a, w);
-            rftfsub(n, a, nc, w + nw);
-        } else if (n == 4) {
-            cftfsub(n, a, w);
-        }
-        xr = a[0] - a[1];
-        a[0] += a[1];
-        for (j = 2; j < n; j += 2) {
-            a[j - 1] = a[j] - a[j + 1];
-            a[j] += a[j + 1];
-        }
-        a[n - 1] = xr;
-    }
-}
-
-
-static void ddst(int n, int isgn, float *a, int *ip, float *w)
-{
-    int j, nw, nc;
-    float xr;
-
-    nw = ip[0];
-    if (n > (nw << 2)) {
-        nw = n >> 2;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > nc) {
-        nc = n;
-        makect(nc, ip, w + nw);
-    }
-    if (isgn < 0) {
-        xr = a[n - 1];
-        for (j = n - 2; j >= 2; j -= 2) {
-            a[j + 1] = -a[j] - a[j - 1];
-            a[j] -= a[j - 1];
-        }
-        a[1] = a[0] + xr;
-        a[0] -= xr;
-        if (n > 4) {
-            rftbsub(n, a, nc, w + nw);
-            bitrv2(n, ip + 2, a);
-            cftbsub(n, a, w);
-        } else if (n == 4) {
-            cftfsub(n, a, w);
-        }
-    }
-    dstsub(n, a, nc, w + nw);
-    if (isgn >= 0) {
-        if (n > 4) {
-            bitrv2(n, ip + 2, a);
-            cftfsub(n, a, w);
-            rftfsub(n, a, nc, w + nw);
-        } else if (n == 4) {
-            cftfsub(n, a, w);
-        }
-        xr = a[0] - a[1];
-        a[0] += a[1];
-        for (j = 2; j < n; j += 2) {
-            a[j - 1] = -a[j] - a[j + 1];
-            a[j] -= a[j + 1];
-        }
-        a[n - 1] = -xr;
-    }
-}
-
-
-static void dfct(int n, float *a, float *t, int *ip, float *w)
-{
-    int j, k, l, m, mh, nw, nc;
-    float xr, xi, yr, yi;
-
-    nw = ip[0];
-    if (n > (nw << 3)) {
-        nw = n >> 3;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > (nc << 1)) {
-        nc = n >> 1;
-        makect(nc, ip, w + nw);
-    }
-    m = n >> 1;
-    yi = a[m];
-    xi = a[0] + a[n];
-    a[0] -= a[n];
-    t[0] = xi - yi;
-    t[m] = xi + yi;
-    if (n > 2) {
-        mh = m >> 1;
-        for (j = 1; j < mh; j++) {
-            k = m - j;
-            xr = a[j] - a[n - j];
-            xi = a[j] + a[n - j];
-            yr = a[k] - a[n - k];
-            yi = a[k] + a[n - k];
-            a[j] = xr;
-            a[k] = yr;
-            t[j] = xi - yi;
-            t[k] = xi + yi;
-        }
-        t[mh] = a[mh] + a[n - mh];
-        a[mh] -= a[n - mh];
-        dctsub(m, a, nc, w + nw);
-        if (m > 4) {
-            bitrv2(m, ip + 2, a);
-            cftfsub(m, a, w);
-            rftfsub(m, a, nc, w + nw);
-        } else if (m == 4) {
-            cftfsub(m, a, w);
-        }
-        a[n - 1] = a[0] - a[1];
-        a[1] = a[0] + a[1];
-        for (j = m - 2; j >= 2; j -= 2) {
-            a[2 * j + 1] = a[j] + a[j + 1];
-            a[2 * j - 1] = a[j] - a[j + 1];
-        }
-        l = 2;
-        m = mh;
-        while (m >= 2) {
-            dctsub(m, t, nc, w + nw);
-            if (m > 4) {
-                bitrv2(m, ip + 2, t);
-                cftfsub(m, t, w);
-                rftfsub(m, t, nc, w + nw);
-            } else if (m == 4) {
-                cftfsub(m, t, w);
-            }
-            a[n - l] = t[0] - t[1];
-            a[l] = t[0] + t[1];
-            k = 0;
-            for (j = 2; j < m; j += 2) {
-                k += l << 2;
-                a[k - l] = t[j] - t[j + 1];
-                a[k + l] = t[j] + t[j + 1];
-            }
-            l <<= 1;
-            mh = m >> 1;
-            for (j = 0; j < mh; j++) {
-                k = m - j;
-                t[j] = t[m + k] - t[m + j];
-                t[k] = t[m + k] + t[m + j];
-            }
-            t[mh] = t[m + mh];
-            m = mh;
-        }
-        a[l] = t[0];
-        a[n] = t[2] - t[1];
-        a[0] = t[2] + t[1];
-    } else {
-        a[1] = a[0];
-        a[2] = t[0];
-        a[0] = t[1];
-    }
-}
-
-static void dfst(int n, float *a, float *t, int *ip, float *w)
-{
-    int j, k, l, m, mh, nw, nc;
-    float xr, xi, yr, yi;
-
-    nw = ip[0];
-    if (n > (nw << 3)) {
-        nw = n >> 3;
-        makewt(nw, ip, w);
-    }
-    nc = ip[1];
-    if (n > (nc << 1)) {
-        nc = n >> 1;
-        makect(nc, ip, w + nw);
-    }
-    if (n > 2) {
-        m = n >> 1;
-        mh = m >> 1;
-        for (j = 1; j < mh; j++) {
-            k = m - j;
-            xr = a[j] + a[n - j];
-            xi = a[j] - a[n - j];
-            yr = a[k] + a[n - k];
-            yi = a[k] - a[n - k];
-            a[j] = xr;
-            a[k] = yr;
-            t[j] = xi + yi;
-            t[k] = xi - yi;
-        }
-        t[0] = a[mh] - a[n - mh];
-        a[mh] += a[n - mh];
-        a[0] = a[m];
-        dstsub(m, a, nc, w + nw);
-        if (m > 4) {
-            bitrv2(m, ip + 2, a);
-            cftfsub(m, a, w);
-            rftfsub(m, a, nc, w + nw);
-        } else if (m == 4) {
-            cftfsub(m, a, w);
-        }
-        a[n - 1] = a[1] - a[0];
-        a[1] = a[0] + a[1];
-        for (j = m - 2; j >= 2; j -= 2) {
-            a[2 * j + 1] = a[j] - a[j + 1];
-            a[2 * j - 1] = -a[j] - a[j + 1];
-        }
-        l = 2;
-        m = mh;
-        while (m >= 2) {
-            dstsub(m, t, nc, w + nw);
-            if (m > 4) {
-                bitrv2(m, ip + 2, t);
-                cftfsub(m, t, w);
-                rftfsub(m, t, nc, w + nw);
-            } else if (m == 4) {
-                cftfsub(m, t, w);
-            }
-            a[n - l] = t[1] - t[0];
-            a[l] = t[0] + t[1];
-            k = 0;
-            for (j = 2; j < m; j += 2) {
-                k += l << 2;
-                a[k - l] = -t[j] - t[j + 1];
-                a[k + l] = t[j] - t[j + 1];
-            }
-            l <<= 1;
-            mh = m >> 1;
-            for (j = 1; j < mh; j++) {
-                k = m - j;
-                t[j] = t[m + k] + t[m + j];
-                t[k] = t[m + k] - t[m + j];
-            }
-            t[0] = t[m + mh];
-            m = mh;
-        }
-        a[l] = t[0];
-    }
-    a[0] = 0;
-}
-#endif  // Not used.
-
-
-/* -------- initializing routines -------- */
-
-
-#include <math.h>
-
-static void makewt(int nw, int *ip, float *w)
-{
-    int j, nwh;
-    float delta, x, y;
-
-    ip[0] = nw;
-    ip[1] = 1;
-    if (nw > 2) {
-        nwh = nw >> 1;
-        delta = (float)atan(1.0f) / nwh;
-        w[0] = 1;
-        w[1] = 0;
-        w[nwh] = (float)cos(delta * nwh);
-        w[nwh + 1] = w[nwh];
-        if (nwh > 2) {
-            for (j = 2; j < nwh; j += 2) {
-                x = (float)cos(delta * j);
-                y = (float)sin(delta * j);
-                w[j] = x;
-                w[j + 1] = y;
-                w[nw - j] = y;
-                w[nw - j + 1] = x;
-            }
-            bitrv2(nw, ip + 2, w);
-        }
-    }
-}
-
-
-static void makect(int nc, int *ip, float *c)
-{
-    int j, nch;
-    float delta;
-
-    ip[1] = nc;
-    if (nc > 1) {
-        nch = nc >> 1;
-        delta = (float)atan(1.0f) / nch;
-        c[0] = (float)cos(delta * nch);
-        c[nch] = 0.5f * c[0];
-        for (j = 1; j < nch; j++) {
-            c[j] = 0.5f * (float)cos(delta * j);
-            c[nc - j] = 0.5f * (float)sin(delta * j);
-        }
-    }
-}
-
-
-/* -------- child routines -------- */
-
-
-static void bitrv2(int n, int *ip, float *a)
-{
-    int j, j1, k, k1, l, m, m2;
-    float xr, xi, yr, yi;
-
-    ip[0] = 0;
-    l = n;
-    m = 1;
-    while ((m << 3) < l) {
-        l >>= 1;
-        for (j = 0; j < m; j++) {
-            ip[m + j] = ip[j] + l;
-        }
-        m <<= 1;
-    }
-    m2 = 2 * m;
-    if ((m << 3) == l) {
-        for (k = 0; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 2 * j + ip[k];
-                k1 = 2 * k + ip[j];
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 += 2 * m2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 -= m2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 += 2 * m2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-            j1 = 2 * k + m2 + ip[k];
-            k1 = j1 + m2;
-            xr = a[j1];
-            xi = a[j1 + 1];
-            yr = a[k1];
-            yi = a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-        }
-    } else {
-        for (k = 1; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 2 * j + ip[k];
-                k1 = 2 * k + ip[j];
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 += m2;
-                xr = a[j1];
-                xi = a[j1 + 1];
-                yr = a[k1];
-                yi = a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-        }
-    }
-}
-
-
-static void bitrv2conj(int n, int *ip, float *a)
-{
-    int j, j1, k, k1, l, m, m2;
-    float xr, xi, yr, yi;
-
-    ip[0] = 0;
-    l = n;
-    m = 1;
-    while ((m << 3) < l) {
-        l >>= 1;
-        for (j = 0; j < m; j++) {
-            ip[m + j] = ip[j] + l;
-        }
-        m <<= 1;
-    }
-    m2 = 2 * m;
-    if ((m << 3) == l) {
-        for (k = 0; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 2 * j + ip[k];
-                k1 = 2 * k + ip[j];
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 += 2 * m2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 -= m2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 += 2 * m2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-            k1 = 2 * k + ip[k];
-            a[k1 + 1] = -a[k1 + 1];
-            j1 = k1 + m2;
-            k1 = j1 + m2;
-            xr = a[j1];
-            xi = -a[j1 + 1];
-            yr = a[k1];
-            yi = -a[k1 + 1];
-            a[j1] = yr;
-            a[j1 + 1] = yi;
-            a[k1] = xr;
-            a[k1 + 1] = xi;
-            k1 += m2;
-            a[k1 + 1] = -a[k1 + 1];
-        }
-    } else {
-        a[1] = -a[1];
-        a[m2 + 1] = -a[m2 + 1];
-        for (k = 1; k < m; k++) {
-            for (j = 0; j < k; j++) {
-                j1 = 2 * j + ip[k];
-                k1 = 2 * k + ip[j];
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-                j1 += m2;
-                k1 += m2;
-                xr = a[j1];
-                xi = -a[j1 + 1];
-                yr = a[k1];
-                yi = -a[k1 + 1];
-                a[j1] = yr;
-                a[j1 + 1] = yi;
-                a[k1] = xr;
-                a[k1 + 1] = xi;
-            }
-            k1 = 2 * k + ip[k];
-            a[k1 + 1] = -a[k1 + 1];
-            a[k1 + m2 + 1] = -a[k1 + m2 + 1];
-        }
-    }
-}
-
-
-static void cftfsub(int n, float *a, float *w)
-{
-    int j, j1, j2, j3, l;
-    float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
-    l = 2;
-    if (n > 8) {
-        cft1st(n, a, w);
-        l = 8;
-        while ((l << 2) < n) {
-            cftmdl(n, l, a, w);
-            l <<= 2;
-        }
-    }
-    if ((l << 2) == n) {
-        for (j = 0; j < l; j += 2) {
-            j1 = j + l;
-            j2 = j1 + l;
-            j3 = j2 + l;
-            x0r = a[j] + a[j1];
-            x0i = a[j + 1] + a[j1 + 1];
-            x1r = a[j] - a[j1];
-            x1i = a[j + 1] - a[j1 + 1];
-            x2r = a[j2] + a[j3];
-            x2i = a[j2 + 1] + a[j3 + 1];
-            x3r = a[j2] - a[j3];
-            x3i = a[j2 + 1] - a[j3 + 1];
-            a[j] = x0r + x2r;
-            a[j + 1] = x0i + x2i;
-            a[j2] = x0r - x2r;
-            a[j2 + 1] = x0i - x2i;
-            a[j1] = x1r - x3i;
-            a[j1 + 1] = x1i + x3r;
-            a[j3] = x1r + x3i;
-            a[j3 + 1] = x1i - x3r;
-        }
-    } else {
-        for (j = 0; j < l; j += 2) {
-            j1 = j + l;
-            x0r = a[j] - a[j1];
-            x0i = a[j + 1] - a[j1 + 1];
-            a[j] += a[j1];
-            a[j + 1] += a[j1 + 1];
-            a[j1] = x0r;
-            a[j1 + 1] = x0i;
-        }
-    }
-}
-
-
-static void cftbsub(int n, float *a, float *w)
-{
-    int j, j1, j2, j3, l;
-    float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
-    l = 2;
-    if (n > 8) {
-        cft1st(n, a, w);
-        l = 8;
-        while ((l << 2) < n) {
-            cftmdl(n, l, a, w);
-            l <<= 2;
-        }
-    }
-    if ((l << 2) == n) {
-        for (j = 0; j < l; j += 2) {
-            j1 = j + l;
-            j2 = j1 + l;
-            j3 = j2 + l;
-            x0r = a[j] + a[j1];
-            x0i = -a[j + 1] - a[j1 + 1];
-            x1r = a[j] - a[j1];
-            x1i = -a[j + 1] + a[j1 + 1];
-            x2r = a[j2] + a[j3];
-            x2i = a[j2 + 1] + a[j3 + 1];
-            x3r = a[j2] - a[j3];
-            x3i = a[j2 + 1] - a[j3 + 1];
-            a[j] = x0r + x2r;
-            a[j + 1] = x0i - x2i;
-            a[j2] = x0r - x2r;
-            a[j2 + 1] = x0i + x2i;
-            a[j1] = x1r - x3i;
-            a[j1 + 1] = x1i - x3r;
-            a[j3] = x1r + x3i;
-            a[j3 + 1] = x1i + x3r;
-        }
-    } else {
-        for (j = 0; j < l; j += 2) {
-            j1 = j + l;
-            x0r = a[j] - a[j1];
-            x0i = -a[j + 1] + a[j1 + 1];
-            a[j] += a[j1];
-            a[j + 1] = -a[j + 1] - a[j1 + 1];
-            a[j1] = x0r;
-            a[j1 + 1] = x0i;
-        }
-    }
-}
-
-
-static void cft1st(int n, float *a, float *w)
-{
-    int j, k1, k2;
-    float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
-    float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
-    x0r = a[0] + a[2];
-    x0i = a[1] + a[3];
-    x1r = a[0] - a[2];
-    x1i = a[1] - a[3];
-    x2r = a[4] + a[6];
-    x2i = a[5] + a[7];
-    x3r = a[4] - a[6];
-    x3i = a[5] - a[7];
-    a[0] = x0r + x2r;
-    a[1] = x0i + x2i;
-    a[4] = x0r - x2r;
-    a[5] = x0i - x2i;
-    a[2] = x1r - x3i;
-    a[3] = x1i + x3r;
-    a[6] = x1r + x3i;
-    a[7] = x1i - x3r;
-    wk1r = w[2];
-    x0r = a[8] + a[10];
-    x0i = a[9] + a[11];
-    x1r = a[8] - a[10];
-    x1i = a[9] - a[11];
-    x2r = a[12] + a[14];
-    x2i = a[13] + a[15];
-    x3r = a[12] - a[14];
-    x3i = a[13] - a[15];
-    a[8] = x0r + x2r;
-    a[9] = x0i + x2i;
-    a[12] = x2i - x0i;
-    a[13] = x0r - x2r;
-    x0r = x1r - x3i;
-    x0i = x1i + x3r;
-    a[10] = wk1r * (x0r - x0i);
-    a[11] = wk1r * (x0r + x0i);
-    x0r = x3i + x1r;
-    x0i = x3r - x1i;
-    a[14] = wk1r * (x0i - x0r);
-    a[15] = wk1r * (x0i + x0r);
-    k1 = 0;
-    for (j = 16; j < n; j += 16) {
-        k1 += 2;
-        k2 = 2 * k1;
-        wk2r = w[k1];
-        wk2i = w[k1 + 1];
-        wk1r = w[k2];
-        wk1i = w[k2 + 1];
-        wk3r = wk1r - 2 * wk2i * wk1i;
-        wk3i = 2 * wk2i * wk1r - wk1i;
-        x0r = a[j] + a[j + 2];
-        x0i = a[j + 1] + a[j + 3];
-        x1r = a[j] - a[j + 2];
-        x1i = a[j + 1] - a[j + 3];
-        x2r = a[j + 4] + a[j + 6];
-        x2i = a[j + 5] + a[j + 7];
-        x3r = a[j + 4] - a[j + 6];
-        x3i = a[j + 5] - a[j + 7];
-        a[j] = x0r + x2r;
-        a[j + 1] = x0i + x2i;
-        x0r -= x2r;
-        x0i -= x2i;
-        a[j + 4] = wk2r * x0r - wk2i * x0i;
-        a[j + 5] = wk2r * x0i + wk2i * x0r;
-        x0r = x1r - x3i;
-        x0i = x1i + x3r;
-        a[j + 2] = wk1r * x0r - wk1i * x0i;
-        a[j + 3] = wk1r * x0i + wk1i * x0r;
-        x0r = x1r + x3i;
-        x0i = x1i - x3r;
-        a[j + 6] = wk3r * x0r - wk3i * x0i;
-        a[j + 7] = wk3r * x0i + wk3i * x0r;
-        wk1r = w[k2 + 2];
-        wk1i = w[k2 + 3];
-        wk3r = wk1r - 2 * wk2r * wk1i;
-        wk3i = 2 * wk2r * wk1r - wk1i;
-        x0r = a[j + 8] + a[j + 10];
-        x0i = a[j + 9] + a[j + 11];
-        x1r = a[j + 8] - a[j + 10];
-        x1i = a[j + 9] - a[j + 11];
-        x2r = a[j + 12] + a[j + 14];
-        x2i = a[j + 13] + a[j + 15];
-        x3r = a[j + 12] - a[j + 14];
-        x3i = a[j + 13] - a[j + 15];
-        a[j + 8] = x0r + x2r;
-        a[j + 9] = x0i + x2i;
-        x0r -= x2r;
-        x0i -= x2i;
-        a[j + 12] = -wk2i * x0r - wk2r * x0i;
-        a[j + 13] = -wk2i * x0i + wk2r * x0r;
-        x0r = x1r - x3i;
-        x0i = x1i + x3r;
-        a[j + 10] = wk1r * x0r - wk1i * x0i;
-        a[j + 11] = wk1r * x0i + wk1i * x0r;
-        x0r = x1r + x3i;
-        x0i = x1i - x3r;
-        a[j + 14] = wk3r * x0r - wk3i * x0i;
-        a[j + 15] = wk3r * x0i + wk3i * x0r;
-    }
-}
-
-
-static void cftmdl(int n, int l, float *a, float *w)
-{
-    int j, j1, j2, j3, k, k1, k2, m, m2;
-    float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
-    float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
-
-    m = l << 2;
-    for (j = 0; j < l; j += 2) {
-        j1 = j + l;
-        j2 = j1 + l;
-        j3 = j2 + l;
-        x0r = a[j] + a[j1];
-        x0i = a[j + 1] + a[j1 + 1];
-        x1r = a[j] - a[j1];
-        x1i = a[j + 1] - a[j1 + 1];
-        x2r = a[j2] + a[j3];
-        x2i = a[j2 + 1] + a[j3 + 1];
-        x3r = a[j2] - a[j3];
-        x3i = a[j2 + 1] - a[j3 + 1];
-        a[j] = x0r + x2r;
-        a[j + 1] = x0i + x2i;
-        a[j2] = x0r - x2r;
-        a[j2 + 1] = x0i - x2i;
-        a[j1] = x1r - x3i;
-        a[j1 + 1] = x1i + x3r;
-        a[j3] = x1r + x3i;
-        a[j3 + 1] = x1i - x3r;
-    }
-    wk1r = w[2];
-    for (j = m; j < l + m; j += 2) {
-        j1 = j + l;
-        j2 = j1 + l;
-        j3 = j2 + l;
-        x0r = a[j] + a[j1];
-        x0i = a[j + 1] + a[j1 + 1];
-        x1r = a[j] - a[j1];
-        x1i = a[j + 1] - a[j1 + 1];
-        x2r = a[j2] + a[j3];
-        x2i = a[j2 + 1] + a[j3 + 1];
-        x3r = a[j2] - a[j3];
-        x3i = a[j2 + 1] - a[j3 + 1];
-        a[j] = x0r + x2r;
-        a[j + 1] = x0i + x2i;
-        a[j2] = x2i - x0i;
-        a[j2 + 1] = x0r - x2r;
-        x0r = x1r - x3i;
-        x0i = x1i + x3r;
-        a[j1] = wk1r * (x0r - x0i);
-        a[j1 + 1] = wk1r * (x0r + x0i);
-        x0r = x3i + x1r;
-        x0i = x3r - x1i;
-        a[j3] = wk1r * (x0i - x0r);
-        a[j3 + 1] = wk1r * (x0i + x0r);
-    }
-    k1 = 0;
-    m2 = 2 * m;
-    for (k = m2; k < n; k += m2) {
-        k1 += 2;
-        k2 = 2 * k1;
-        wk2r = w[k1];
-        wk2i = w[k1 + 1];
-        wk1r = w[k2];
-        wk1i = w[k2 + 1];
-        wk3r = wk1r - 2 * wk2i * wk1i;
-        wk3i = 2 * wk2i * wk1r - wk1i;
-        for (j = k; j < l + k; j += 2) {
-            j1 = j + l;
-            j2 = j1 + l;
-            j3 = j2 + l;
-            x0r = a[j] + a[j1];
-            x0i = a[j + 1] + a[j1 + 1];
-            x1r = a[j] - a[j1];
-            x1i = a[j + 1] - a[j1 + 1];
-            x2r = a[j2] + a[j3];
-            x2i = a[j2 + 1] + a[j3 + 1];
-            x3r = a[j2] - a[j3];
-            x3i = a[j2 + 1] - a[j3 + 1];
-            a[j] = x0r + x2r;
-            a[j + 1] = x0i + x2i;
-            x0r -= x2r;
-            x0i -= x2i;
-            a[j2] = wk2r * x0r - wk2i * x0i;
-            a[j2 + 1] = wk2r * x0i + wk2i * x0r;
-            x0r = x1r - x3i;
-            x0i = x1i + x3r;
-            a[j1] = wk1r * x0r - wk1i * x0i;
-            a[j1 + 1] = wk1r * x0i + wk1i * x0r;
-            x0r = x1r + x3i;
-            x0i = x1i - x3r;
-            a[j3] = wk3r * x0r - wk3i * x0i;
-            a[j3 + 1] = wk3r * x0i + wk3i * x0r;
-        }
-        wk1r = w[k2 + 2];
-        wk1i = w[k2 + 3];
-        wk3r = wk1r - 2 * wk2r * wk1i;
-        wk3i = 2 * wk2r * wk1r - wk1i;
-        for (j = k + m; j < l + (k + m); j += 2) {
-            j1 = j + l;
-            j2 = j1 + l;
-            j3 = j2 + l;
-            x0r = a[j] + a[j1];
-            x0i = a[j + 1] + a[j1 + 1];
-            x1r = a[j] - a[j1];
-            x1i = a[j + 1] - a[j1 + 1];
-            x2r = a[j2] + a[j3];
-            x2i = a[j2 + 1] + a[j3 + 1];
-            x3r = a[j2] - a[j3];
-            x3i = a[j2 + 1] - a[j3 + 1];
-            a[j] = x0r + x2r;
-            a[j + 1] = x0i + x2i;
-            x0r -= x2r;
-            x0i -= x2i;
-            a[j2] = -wk2i * x0r - wk2r * x0i;
-            a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
-            x0r = x1r - x3i;
-            x0i = x1i + x3r;
-            a[j1] = wk1r * x0r - wk1i * x0i;
-            a[j1 + 1] = wk1r * x0i + wk1i * x0r;
-            x0r = x1r + x3i;
-            x0i = x1i - x3r;
-            a[j3] = wk3r * x0r - wk3i * x0i;
-            a[j3 + 1] = wk3r * x0i + wk3i * x0r;
-        }
-    }
-}
-
-
-static void rftfsub(int n, float *a, int nc, float *c)
-{
-    int j, k, kk, ks, m;
-    float wkr, wki, xr, xi, yr, yi;
-
-    m = n >> 1;
-    ks = 2 * nc / m;
-    kk = 0;
-    for (j = 2; j < m; j += 2) {
-        k = n - j;
-        kk += ks;
-        wkr = 0.5f - c[nc - kk];
-        wki = c[kk];
-        xr = a[j] - a[k];
-        xi = a[j + 1] + a[k + 1];
-        yr = wkr * xr - wki * xi;
-        yi = wkr * xi + wki * xr;
-        a[j] -= yr;
-        a[j + 1] -= yi;
-        a[k] += yr;
-        a[k + 1] -= yi;
-    }
-}
-
-
-static void rftbsub(int n, float *a, int nc, float *c)
-{
-    int j, k, kk, ks, m;
-    float wkr, wki, xr, xi, yr, yi;
-
-    a[1] = -a[1];
-    m = n >> 1;
-    ks = 2 * nc / m;
-    kk = 0;
-    for (j = 2; j < m; j += 2) {
-        k = n - j;
-        kk += ks;
-        wkr = 0.5f - c[nc - kk];
-        wki = c[kk];
-        xr = a[j] - a[k];
-        xi = a[j + 1] + a[k + 1];
-        yr = wkr * xr + wki * xi;
-        yi = wkr * xi - wki * xr;
-        a[j] -= yr;
-        a[j + 1] = yi - a[j + 1];
-        a[k] += yr;
-        a[k + 1] = yi - a[k + 1];
-    }
-    a[m + 1] = -a[m + 1];
-}
-
-#if 0  // Not used.
-static void dctsub(int n, float *a, int nc, float *c)
-{
-    int j, k, kk, ks, m;
-    float wkr, wki, xr;
-
-    m = n >> 1;
-    ks = nc / n;
-    kk = 0;
-    for (j = 1; j < m; j++) {
-        k = n - j;
-        kk += ks;
-        wkr = c[kk] - c[nc - kk];
-        wki = c[kk] + c[nc - kk];
-        xr = wki * a[j] - wkr * a[k];
-        a[j] = wkr * a[j] + wki * a[k];
-        a[k] = xr;
-    }
-    a[m] *= c[0];
-}
-
-
-static void dstsub(int n, float *a, int nc, float *c)
-{
-    int j, k, kk, ks, m;
-    float wkr, wki, xr;
-
-    m = n >> 1;
-    ks = nc / n;
-    kk = 0;
-    for (j = 1; j < m; j++) {
-        k = n - j;
-        kk += ks;
-        wkr = c[kk] - c[nc - kk];
-        wki = c[kk] + c[nc - kk];
-        xr = wki * a[k] - wkr * a[j];
-        a[k] = wkr * a[k] + wki * a[j];
-        a[j] = xr;
-    }
-    a[m] *= c[0];
-}
-#endif  // Not used.
diff --git a/webrtc/modules/audio_processing/utility/fft4g.h b/webrtc/modules/audio_processing/utility/fft4g.h
deleted file mode 100644
index 14a52a106a..0000000000
--- a/webrtc/modules/audio_processing/utility/fft4g.h
+++ /dev/null
@@ -1,17 +0,0 @@
-/*
- *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
- *
- *  Use of this source code is governed by a BSD-style license
- *  that can be found in the LICENSE file in the root of the source
- *  tree. An additional intellectual property rights grant can be found
- *  in the file PATENTS.  All contributing project authors may
- *  be found in the AUTHORS file in the root of the source tree.
- */
-
-#ifndef WEBRTC_MODULES_AUDIO_PROCESSING_UTILITY_FFT4G_H_
-#define WEBRTC_MODULES_AUDIO_PROCESSING_UTILITY_FFT4G_H_
-
-void WebRtc_rdft(int, int, float *, int *, float *);
-void WebRtc_cdft(int, int, float *, int *, float *);
-
-#endif