ref: e29dcf169ba38a15396dafcd40e6d405e5d4abd7
parent: 3e7c1e463cc93eb7b8f3aae049a20c6dc1ac3e0b
author: robs <robs>
date: Sat Sep 13 07:59:38 EDT 2008
fft tidy-up
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -28,16 +28,15 @@
# Format with: !xargs echo|tr ' ' '\n'|sort|column|expand|sed 's/^/ /'
set(effects_srcs
- biquad echo noiseprof remix swap
- biquads echos noisered repeat synth
- chorus fade normalise resample tempo
- compand FFT output reverb tremolo
- compandt fft4g pad reverse trim
- contrast filter pan silence vol
- dcshift flanger phaser skeleff
- delay input pitch speed
- dither mcompand polyphas splice
- earwax mixer rate stat
+ biquad earwax mixer polyphas speed
+ biquads echo noiseprof rate splice
+ chorus echos noisered remix stat
+ compand fade normalise repeat swap
+ compandt fft4g output resample synth
+ contrast filter pad reverb tempo
+ dcshift flanger pan reverse tremolo
+ delay input phaser silence trim
+ dither mcompand pitch skeleff vol
)
set(formats_srcs
8svx dat ima-fmt s3-fmt u4-fmt
--- a/src/FFT.c
+++ /dev/null
@@ -1,380 +1,0 @@
-/* libSoX FFT funtions copyright Ian Turner and others.
- *
- * This library is free software; you can redistribute it and/or modify it
- * under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or (at
- * your option) any later version.
- *
- * This library is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
- * General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this library; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
- *
- *
- * Based on FFT.cpp from Audacity, with the following permission from
- * its author, Dominic Mazzoni (in particular, relicensing the code
- * for use in SoX):
- *
- * I hereby license you under the LGPL all of the code in FFT.cpp
- * from any version of Audacity, with the exception of the windowing
- * function code, as I wrote the rest of the line [sic] that appears
- * in any version of Audacity (or they are derived from Don Cross or
- * NR, which is okay).
- *
- * -- Dominic Mazzoni <dominic@audacityteam.org>, 18th November 2006
- *
- * As regards the windowing function, WindowFunc, Dominic granted a
- * license to it too, writing on the same day:
- *
- * OK, we're good. That's the original version that I wrote, before
- * others contributed.
- *
- * Some of this code was based on a free implementation of an FFT
- * by Don Cross, available on the web at:
- *
- * http://www.intersrv.com/~dcross/fft.html [no longer, it seems]
- *
- * The basic algorithm for his code was based on Numerical Recipes
- * in Fortran.
- */
-
-#include "sox_i.h"
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <assert.h>
-
-#include "FFT.h"
-
-unsigned **gFFTBitTable = NULL;
-const unsigned MaxFastBits = 16;
-
-static int IsPowerOfTwo(unsigned x)
-{- if (x < 2)
- return 0;
-
- return !(x & (x-1)); /* Thanks to 'byang' for this cute trick! */
-}
-
-static int NumberOfBitsNeeded(unsigned PowerOfTwo)
-{- int i;
-
- if (PowerOfTwo < 2) {- sox_fail("Error: FFT called with size %d", PowerOfTwo);- exit(2);
- }
-
- for (i = 0;; i++)
- if (PowerOfTwo & (1 << i))
- return i;
-}
-
-static unsigned ReverseBits(unsigned index, unsigned NumBits)
-{- unsigned i, rev;
-
- for (i = rev = 0; i < NumBits; i++) {- rev = (rev << 1) | (index & 1);
- index >>= 1;
- }
-
- return rev;
-}
-
-
-/* This function permanently allocates about 250Kb (actually (2**16)-2 ints). */
-static void InitFFT(void)
-{- unsigned len, b;
-
- gFFTBitTable = lsx_calloc(MaxFastBits, sizeof(*gFFTBitTable));
-
- for (b = 1, len = 2; b <= MaxFastBits; b++) {- unsigned i;
-
- gFFTBitTable[b - 1] = lsx_calloc(len, sizeof(**gFFTBitTable));
- for (i = 0; i < len; i++)
- gFFTBitTable[b - 1][i] = ReverseBits(i, b);
-
- len <<= 1;
- }
-}
-
-#define FastReverseBits(i, NumBits) \
- (NumBits <= MaxFastBits) ? gFFTBitTable[NumBits - 1][i] : ReverseBits(i, NumBits)
-
-/*
- * Complex Fast Fourier Transform
- */
-void FFT(unsigned NumSamples,
- int InverseTransform,
- const float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
-{- unsigned i, BlockSize, NumBits; /* Number of bits needed to store indices */
- int j, k, n;
- int BlockEnd;
-
- double angle_numerator = 2.0 * M_PI;
- double tr, ti; /* temp real, temp imaginary */
-
- if (!IsPowerOfTwo(NumSamples)) {- sox_fail("%d is not a power of two", NumSamples);- exit(2);
- }
-
- if (!gFFTBitTable)
- InitFFT();
-
- if (InverseTransform)
- angle_numerator = -angle_numerator;
-
- NumBits = NumberOfBitsNeeded(NumSamples);
-
- /*
- ** Do simultaneous data copy and bit-reversal ordering into outputs...
- */
-
- for (i = 0; i < NumSamples; i++) {- j = FastReverseBits(i, NumBits);
- RealOut[j] = RealIn[i];
- ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
- }
-
- /*
- ** Do the FFT itself...
- */
-
- BlockEnd = 1;
- for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {-
- double delta_angle = angle_numerator / (double) BlockSize;
-
- double sm2 = sin(-2 * delta_angle);
- double sm1 = sin(-delta_angle);
- double cm2 = cos(-2 * delta_angle);
- double cm1 = cos(-delta_angle);
- double w = 2 * cm1;
- double ar0, ar1, ar2, ai0, ai1, ai2;
-
- for (i = 0; i < NumSamples; i += BlockSize) {- ar2 = cm2;
- ar1 = cm1;
-
- ai2 = sm2;
- ai1 = sm1;
-
- for (j = i, n = 0; n < BlockEnd; j++, n++) {- ar0 = w * ar1 - ar2;
- ar2 = ar1;
- ar1 = ar0;
-
- ai0 = w * ai1 - ai2;
- ai2 = ai1;
- ai1 = ai0;
-
- k = j + BlockEnd;
- tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
- ti = ar0 * ImagOut[k] + ai0 * RealOut[k];
-
- RealOut[k] = RealOut[j] - tr;
- ImagOut[k] = ImagOut[j] - ti;
-
- RealOut[j] += tr;
- ImagOut[j] += ti;
- }
- }
-
- BlockEnd = BlockSize;
- }
-
- /*
- ** Need to normalize if inverse transform...
- */
-
- if (InverseTransform) {- float denom = (float) NumSamples;
-
- for (i = 0; i < NumSamples; i++) {- RealOut[i] /= denom;
- ImagOut[i] /= denom;
- }
- }
-}
-
-/*
- * Real Fast Fourier Transform
- *
- * This function was based on the code in Numerical Recipes for C.
- * In Num. Rec., the inner loop is based on a single 1-based array
- * of interleaved real and imaginary numbers. Because we have two
- * separate zero-based arrays, our indices are quite different.
- * Here is the correspondence between Num. Rec. indices and our indices:
- *
- * i1 <-> real[i]
- * i2 <-> imag[i]
- * i3 <-> real[n/2-i]
- * i4 <-> imag[n/2-i]
- */
-
-void RealFFT(unsigned NumSamples, const float *RealIn, float *RealOut, float *ImagOut)
-{- unsigned i, i3, Half = NumSamples / 2;
- float theta = M_PI / Half;
- float wtemp = (float) sin(0.5 * theta);
- float wpr = -2.0 * wtemp * wtemp;
- float wpi = (float) sin(theta);
- float wr = 1.0 + wpr;
- float wi = wpi;
- float h1r, h1i, h2r, h2i;
- float *tmpReal, *tmpImag;
-
- tmpReal = lsx_calloc(NumSamples, sizeof(float));
- tmpImag = tmpReal + Half;
-
- for (i = 0; i < Half; i++) {- tmpReal[i] = RealIn[2 * i];
- tmpImag[i] = RealIn[2 * i + 1];
- }
-
- FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
-
- for (i = 1; i < Half / 2; i++) {- i3 = Half - i;
-
- h1r = 0.5 * (RealOut[i] + RealOut[i3]);
- h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
- h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
- h2i = -0.5 * (RealOut[i] - RealOut[i3]);
-
- RealOut[i] = h1r + wr * h2r - wi * h2i;
- ImagOut[i] = h1i + wr * h2i + wi * h2r;
- RealOut[i3] = h1r - wr * h2r + wi * h2i;
- ImagOut[i3] = -h1i + wr * h2i + wi * h2r;
-
- wtemp = wr;
- wr = wr * wpr - wi * wpi + wr;
- wi = wi * wpr + wtemp * wpi + wi;
- }
-
- h1r = RealOut[0];
- RealOut[0] += ImagOut[0];
- ImagOut[0] = h1r - ImagOut[0];
-
- free(tmpReal);
-}
-
-/*
- * PowerSpectrum
- *
- * This function computes the same as RealFFT, above, but
- * adds the squares of the real and imaginary part of each
- * coefficient, extracting the power and throwing away the
- * phase.
- *
- * For speed, it does not call RealFFT, but duplicates some
- * of its code.
- */
-
-void PowerSpectrum(unsigned NumSamples, const float *In, float *Out)
-{- unsigned Half, i, i3;
- float theta, wtemp, wpr, wpi, wr, wi;
- float h1r, h1i, h2r, h2i, rt, it;
- float *tmpReal;
- float *tmpImag;
- float *RealOut;
- float *ImagOut;
-
- Half = NumSamples / 2;
-
- theta = M_PI / Half;
-
- tmpReal = lsx_calloc(Half * 4, sizeof(float));
- tmpImag = tmpReal + Half;
- RealOut = tmpImag + Half;
- ImagOut = RealOut + Half;
-
- for (i = 0; i < Half; i++) {- tmpReal[i] = In[2 * i];
- tmpImag[i] = In[2 * i + 1];
- }
-
- FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
-
- wtemp = (float) sin(0.5 * theta);
-
- wpr = -2.0 * wtemp * wtemp;
- wpi = (float) sin(theta);
- wr = 1.0 + wpr;
- wi = wpi;
-
- for (i = 1; i < Half / 2; i++) {-
- i3 = Half - i;
-
- h1r = 0.5 * (RealOut[i] + RealOut[i3]);
- h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
- h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
- h2i = -0.5 * (RealOut[i] - RealOut[i3]);
-
- rt = h1r + wr * h2r - wi * h2i;
- it = h1i + wr * h2i + wi * h2r;
-
- Out[i] = rt * rt + it * it;
-
- rt = h1r - wr * h2r + wi * h2i;
- it = -h1i + wr * h2i + wi * h2r;
-
- Out[i3] = rt * rt + it * it;
-
- wr = (wtemp = wr) * wpr - wi * wpi + wr;
- wi = wi * wpr + wtemp * wpi + wi;
- }
-
- rt = (h1r = RealOut[0]) + ImagOut[0];
- it = h1r - ImagOut[0];
- Out[0] = rt * rt + it * it;
-
- rt = RealOut[Half / 2];
- it = ImagOut[Half / 2];
- Out[Half / 2] = rt * rt + it * it;
-
- free(tmpReal);
-}
-
-/*
- * Windowing Functions
- */
-
-void WindowFunc(windowfunc_t whichFunction, int NumSamples, float *in)
-{- int i;
-
- switch (whichFunction) {- case BARTLETT:
- for (i = 0; i < NumSamples / 2; i++) {- in[i] *= (i / (float) (NumSamples / 2));
- in[i + (NumSamples / 2)] *=
- (1.0 - (i / (float) (NumSamples / 2)));
- }
- break;
-
- case HAMMING:
- for (i = 0; i < NumSamples; i++)
- in[i] *= 0.54 - 0.46 * cos(2 * M_PI * i / (NumSamples - 1));
- break;
-
- case HANNING:
- for (i = 0; i < NumSamples; i++)
- in[i] *= 0.50 - 0.50 * cos(2 * M_PI * i / (NumSamples - 1));
- break;
- case RECTANGULAR:
- break;
- }
-}
--- a/src/FFT.h
+++ /dev/null
@@ -1,92 +1,0 @@
-/* libSoX FFT funtions copyright Ian Turner and others.
- *
- * This library is free software; you can redistribute it and/or modify it
- * under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or (at
- * your option) any later version.
- *
- * This library is distributed in the hope that it will be useful, but
- * WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
- * General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this library; if not, write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
- *
- *
- * Based on FFT.h from Audacity, with the following permission from
- * its author, Dominic Mazzoni (in particular, relicensing the code
- * for use in SoX):
- *
- * I hereby license you under the LGPL all of the code in FFT.cpp
- * from any version of Audacity, with the exception of the windowing
- * function code, as I wrote the rest of the line [sic] that appears
- * in any version of Audacity (or they are derived from Don Cross or
- * NR, which is okay).
- *
- * -- Dominic Mazzoni <dominic@audacityteam.org>, 18th November 2006
- *
- * As regards the windowing function, WindowFunc, Dominic granted a
- * license to it too, writing on the same day:
- *
- * OK, we're good. That's the original version that I wrote, before
- * others contributed.
- *
- * Some of this code was based on a free implementation of an FFT
- * by Don Cross, available on the web at:
- *
- * http://www.intersrv.com/~dcross/fft.html [no longer, it seems]
- *
- * The basic algorithm for his code was based on Numerical Recipes
- * in Fortran.
- */
-
-/* aliases */
-#define gFFTBitTable lsx_gFFTBitTable
-#define MaxFastBits lsx_MaxFastBits
-#define FFT lsx_FFT
-#define RealFFT lsx_RealFFT
-#define PowerSpectrum lsx_PowerSpectrum
-#define WindowFunc lsx_WindowFunc
-
-/*
- * This is the function you will use the most often.
- * Given an array of floats, this will compute the power
- * spectrum by doing a Real FFT and then computing the
- * sum of the squares of the real and imaginary parts.
- * Note that the output array is half the length of the
- * input array, and that NumSamples must be a power of two.
- */
-
-void PowerSpectrum(unsigned NumSamples, const float *In, float *Out);
-
-/*
- * Computes an FFT when the input data is real but you still
- * want complex data as output. The output arrays are half
- * the length of the input, and NumSamples must be a power of
- * two.
- */
-
-void RealFFT(unsigned NumSamples,
- const float *RealIn, float *RealOut, float *ImagOut);
-
-/*
- * Computes a FFT of complex input and returns complex output.
- * Currently this is the only function here that supports the
- * inverse transform as well.
- */
-
-void FFT(unsigned NumSamples,
- int InverseTransform,
- const float *RealIn, float *ImagIn, float *RealOut, float *ImagOut);
-
-/*
- * Applies a windowing function to the data in place
- */
-typedef enum {RECTANGULAR, /* no window */- BARTLETT, /* triangular */
- HAMMING,
- HANNING} windowfunc_t;
-
-void WindowFunc(windowfunc_t whichFunction, int NumSamples, float *data);
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -251,8 +251,8 @@
libsfx_la_SOURCES = \
band.h biquad.c biquad.h biquads.c chorus.c compand.c compandt.c \
compandt.h contrast.c dcshift.c delay.c dither.c earwax.c echo.c \
- echos.c effects.c effects.h effects_i.c fade.c fft4g.c fft4g.h FFT.c \
- FFT.h fifo.h filter.c flanger.c input.c ladspa.c mcompand.c \
+ echos.c effects.c effects.h effects_i.c fade.c fft4g.c fft4g.h \
+ fifo.h filter.c flanger.c input.c ladspa.c mcompand.c \
mixer.c noiseprof.c noisered.c noisered.h normalise.c output.c pad.c \
pan.c phaser.c pitch.c polyphas.c rabbit.c rate.c \
rate_filters.h rate_half_fir.h rate_poly_fir0.h rate_poly_fir.h \
--- a/src/effects.c
+++ b/src/effects.c
@@ -83,7 +83,7 @@
int sox_effect_options(sox_effect_t *effp, int argc, char * const argv[])
{- return effp->handler.getopts(effp, argc, argv);
+ return effp->handler.getopts(effp, argc, (char * *)argv);
} /* sox_effect_options */
/* Effects chain: */
--- a/src/effects_i.c
+++ b/src/effects_i.c
@@ -248,7 +248,7 @@
return result < 0 ? -1 : result;
}
-double bessel_I_0(double x)
+double lsx_bessel_I_0(double x)
{double term = 1, sum = 1, last_sum, x2 = x / 2;
int i = 1;
@@ -258,3 +258,126 @@
} while (sum != last_sum);
return sum;
}
+
+#include "fft4g.h"
+int * lsx_fft_br;
+double * lsx_fft_sc;
+
+static void update_fft_cache(int len)
+{+ static int n;
+
+ if (!len) {+ free(lsx_fft_br);
+ free(lsx_fft_sc);
+ lsx_fft_sc = NULL;
+ lsx_fft_br = NULL;
+ n = 0;
+ return;
+ }
+ if (len > n) {+ int old_n = n;
+ n = len;
+ lsx_fft_br = lsx_realloc(lsx_fft_br, dft_br_len(n) * sizeof(*lsx_fft_br));
+ lsx_fft_sc = lsx_realloc(lsx_fft_sc, dft_sc_len(n) * sizeof(*lsx_fft_sc));
+ if (!old_n)
+ lsx_fft_br[0] = 0;
+ }
+}
+
+static sox_bool is_power_of_2(int x)
+{+ return !(x < 2 || (x & (x - 1)));
+}
+
+void lsx_safe_rdft(int len, int type, double * d)
+{+ assert(is_power_of_2(len));
+ update_fft_cache(len);
+ lsx_rdft(len, type, d, lsx_fft_br, lsx_fft_sc);
+}
+
+void lsx_safe_cdft(int len, int type, double * d)
+{+ assert(is_power_of_2(len));
+ update_fft_cache(len);
+ lsx_cdft(len, type, d, lsx_fft_br, lsx_fft_sc);
+}
+
+void lsx_power_spectrum(int n, double const * in, double * out)
+{+ int i;
+ double * work = lsx_memdup(in, n * sizeof(*work));
+ lsx_safe_rdft(n, 1, work);
+ out[0] = sqr(work[0]);
+ for (i = 2; i < n; i += 2)
+ out[i >> 1] = sqr(work[i]) + sqr(work[i + 1]);
+ out[i >> 1] = sqr(work[1]);
+ free(work);
+}
+
+void lsx_power_spectrum_f(int n, float const * in, float * out)
+{+ int i;
+ double * work = lsx_malloc(n * sizeof(*work));
+ for (i = 0; i< n; ++i) work[i] = in[i];
+ lsx_safe_rdft(n, 1, work);
+ out[0] = sqr(work[0]);
+ for (i = 2; i < n; i += 2)
+ out[i >> 1] = sqr(work[i]) + sqr(work[i + 1]);
+ out[i >> 1] = sqr(work[1]);
+ free(work);
+}
+
+void lsx_apply_hann_f(float h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= .5 - .5 * cos(x);
+ }
+}
+
+void lsx_apply_hann(double h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= .5 - .5 * cos(x);
+ }
+}
+
+void lsx_apply_hamming(double h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ double x = 2 * M_PI * i / m;
+ h[i] *= .53836 - .46164 * cos(x);
+ }
+}
+
+void lsx_apply_bartlett(double h[], const int num_points)
+{+ int i, m = num_points - 1;
+ for (i = 0; i < num_points; ++i) {+ h[i] *= 2. / m * (m / 2. - fabs(i - m / 2.));
+ }
+}
+
+double lsx_kaiser_beta(double att)
+{+ if (att > 100 ) return .1117 * att - 1.11;
+ if (att > 50 ) return .1102 * (att - 8.7);
+ if (att > 20.96) return .58417 * pow(att -20.96, .4) + .07886 * (att - 20.96);
+ return 0;
+}
+
+void lsx_apply_kaiser(double h[], const int num_points, double beta)
+{+ int i, m = num_points - 1;
+ for (i = 0; i <= m; ++i) {+ double x = 2. * i / m - 1;
+ h[i] *= lsx_bessel_I_0(beta * sqrt(1 - x * x)) / lsx_bessel_I_0(beta);
+ }
+}
+
--- a/src/fft4g.c
+++ b/src/fft4g.c
@@ -291,12 +291,19 @@
#define cos cosf
#define atan atanf
- #define cdft cdft_f
- #define rdft rdft_f
- #define ddct ddct_f
- #define ddst ddst_f
- #define dfct dfct_f
- #define dfst dfst_f
+ #define cdft lsx_cdft_f
+ #define rdft lsx_rdft_f
+ #define ddct lsx_ddct_f
+ #define ddst lsx_ddst_f
+ #define dfct lsx_dfct_f
+ #define dfst lsx_dfst_f
+#else
+ #define cdft lsx_cdft
+ #define rdft lsx_rdft
+ #define ddct lsx_ddct
+ #define ddst lsx_ddst
+ #define dfct lsx_dfct
+ #define dfst lsx_dfst
#endif
static void bitrv2conj(int n, int *ip, double *a);
--- a/src/fft4g.h
+++ b/src/fft4g.h
@@ -13,19 +13,19 @@
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
-void cdft(int, int, double *, int *, double *);
-void rdft(int, int, double *, int *, double *);
-void ddct(int, int, double *, int *, double *);
-void ddst(int, int, double *, int *, double *);
-void dfct(int, double *, double *, int *, double *);
-void dfst(int, double *, double *, int *, double *);
+void lsx_cdft(int, int, double *, int *, double *);
+void lsx_rdft(int, int, double *, int *, double *);
+void lsx_ddct(int, int, double *, int *, double *);
+void lsx_ddst(int, int, double *, int *, double *);
+void lsx_dfct(int, double *, double *, int *, double *);
+void lsx_dfst(int, double *, double *, int *, double *);
-void cdft_f(int, int, float *, int *, float *);
-void rdft_f(int, int, float *, int *, float *);
-void ddct_f(int, int, float *, int *, float *);
-void ddst_f(int, int, float *, int *, float *);
-void dfct_f(int, float *, float *, int *, float *);
-void dfst_f(int, float *, float *, int *, float *);
+void lsx_cdft_f(int, int, float *, int *, float *);
+void lsx_rdft_f(int, int, float *, int *, float *);
+void lsx_ddct_f(int, int, float *, int *, float *);
+void lsx_ddst_f(int, int, float *, int *, float *);
+void lsx_dfct_f(int, float *, float *, int *, float *);
+void lsx_dfst_f(int, float *, float *, int *, float *);
#define dft_br_len(l) (2 + (1 << (int)(log(l / 2 + .5) / log(2.)) / 2))
#define dft_sc_len(l) (l / 2)
--- a/src/noiseprof.c
+++ b/src/noiseprof.c
@@ -95,7 +95,7 @@
float *out = lsx_calloc(FREQCOUNT, sizeof(float));
int i;
- PowerSpectrum(WINDOWSIZE, chan->window, out);
+ lsx_power_spectrum_f(WINDOWSIZE, chan->window, out);
for (i = 0; i < FREQCOUNT; i ++) { if (out[i] > 0) {--- a/src/noisered.c
+++ b/src/noisered.c
@@ -31,6 +31,28 @@
size_t bufdata;
} priv_t;
+static void FFT(unsigned NumSamples,
+ int InverseTransform,
+ const float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
+{+ unsigned i;
+ double * work = malloc(2 * NumSamples * sizeof(*work));
+ for (i = 0; i < 2 * NumSamples; i += 2) {+ work[i] = RealIn[i >> 1];
+ work[i + 1] = ImagIn? ImagIn[i >> 1] : 0;
+ }
+ lsx_safe_cdft(2 * (int)NumSamples, InverseTransform? -1 : 1, work);
+ if (InverseTransform) for (i = 0; i < 2 * NumSamples; i += 2) {+ RealOut[i >> 1] = work[i] / NumSamples;
+ ImagOut[i >> 1] = work[i + 1] / NumSamples;
+ }
+ else for (i = 0; i < 2 * NumSamples; i += 2) {+ RealOut[i >> 1] = work[i];
+ ImagOut[i >> 1] = work[i + 1];
+ }
+ free(work);
+}
+
/*
* Get the options. Default file is stdin (if the audio
* input file isn't coming from there, of course!)
@@ -146,8 +168,8 @@
FFT(WINDOWSIZE, 0, inr, NULL, outr, outi);
memcpy(inr, window, WINDOWSIZE*sizeof(float));
- WindowFunc(HANNING, WINDOWSIZE, inr);
- PowerSpectrum(WINDOWSIZE, inr, power);
+ lsx_apply_hann_f(inr, WINDOWSIZE);
+ lsx_power_spectrum_f(WINDOWSIZE, inr, power);
for (i = 0; i < FREQCOUNT; i ++) {float smooth;
@@ -189,7 +211,7 @@
}
FFT(WINDOWSIZE, 1, outr, outi, inr, ini);
- WindowFunc(HANNING, WINDOWSIZE, inr);
+ lsx_apply_hann_f(inr, WINDOWSIZE);
memcpy(window, inr, WINDOWSIZE*sizeof(float));
--- a/src/noisered.h
+++ b/src/noisered.h
@@ -19,8 +19,6 @@
*/
#include "sox_i.h"
-#include "FFT.h"
-
#include <math.h>
#define WINDOWSIZE 2048
--- a/src/rate.c
+++ b/src/rate.c
@@ -143,7 +143,7 @@
fifo_trim_by(output_fifo, (f->dft_length + overlap) >> 1);
memcpy(output, input, f->dft_length * sizeof(*output));
- rdft(f->dft_length, 1, output, s->bit_rev_table, s->sin_cos_table);
+ lsx_rdft(f->dft_length, 1, output, s->bit_rev_table, s->sin_cos_table);
output[0] *= f->coefs[0];
output[1] *= f->coefs[1];
for (i = 2; i < f->dft_length; i += 2) {@@ -151,7 +151,7 @@
output[i ] = f->coefs[i ] * tmp - f->coefs[i+1] * output[i+1];
output[i+1] = f->coefs[i+1] * tmp + f->coefs[i ] * output[i+1];
}
- rdft(f->dft_length, -1, output, s->bit_rev_table, s->sin_cos_table);
+ lsx_rdft(f->dft_length, -1, output, s->bit_rev_table, s->sin_cos_table);
for (j = 1, i = 2; i < f->dft_length - overlap; ++j, i += 2)
output[j] = output[i];
@@ -176,7 +176,7 @@
for (j = i = 0; i < f->dft_length; ++j, i += 2)
output[i] = input[j], output[i+1] = 0;
- rdft(f->dft_length, 1, output, s->bit_rev_table, s->sin_cos_table);
+ lsx_rdft(f->dft_length, 1, output, s->bit_rev_table, s->sin_cos_table);
output[0] *= f->coefs[0];
output[1] *= f->coefs[1];
for (i = 2; i < f->dft_length; i += 2) {@@ -184,7 +184,7 @@
output[i ] = f->coefs[i ] * tmp - f->coefs[i+1] * output[i+1];
output[i+1] = f->coefs[i+1] * tmp + f->coefs[i ] * output[i+1];
}
- rdft(f->dft_length, -1, output, s->bit_rev_table, s->sin_cos_table);
+ lsx_rdft(f->dft_length, -1, output, s->bit_rev_table, s->sin_cos_table);
}
}
@@ -196,7 +196,7 @@
for (i = 0; i <= m / 2; ++i) {double x = M_PI * (i - .5 * m), y = 2. * i / m - 1;
h[i] = x? sin(Fc * x) / x : Fc;
- sum += h[i] *= bessel_I_0(beta * sqrt(1 - y * y));
+ sum += h[i] *= lsx_bessel_I_0(beta * sqrt(1 - y * y));
if (m - i != i)
sum += h[m - i] = h[i];
}
@@ -254,19 +254,19 @@
work = calloc(work_len, sizeof(*work));
for (i = 0; i < *len; ++i) work[i] = (*h)[i];
- rdft(work_len, 1, work, p->bit_rev_table, p->sin_cos_table); /* Cepstrum: */
+ lsx_rdft(work_len, 1, work, p->bit_rev_table, p->sin_cos_table); /* Cepstrum: */
work[0] = log(fabs(work[0])), work[1] = log(fabs(work[1]));
for (i = 2; i < work_len; i += 2) {work[i] = log(sqrt(sqr(work[i]) + sqr(work[i + 1])));
work[i + 1] = 0;
}
- rdft(work_len, -1, work, p->bit_rev_table, p->sin_cos_table);
+ lsx_rdft(work_len, -1, work, p->bit_rev_table, p->sin_cos_table);
for (i = 0; i < work_len; ++i) work[i] *= 2. / work_len;
for (i = 1; i < work_len / 2; ++i) { /* Window to reject acausal components */work[i] *= 2;
work[i + work_len / 2] = 0;
}
- rdft(work_len, 1, work, p->bit_rev_table, p->sin_cos_table);
+ lsx_rdft(work_len, 1, work, p->bit_rev_table, p->sin_cos_table);
/* Some DFTs require phase unwrapping here, but rdft seems not to. */
for (i = 2; i < work_len; i += 2) /* Interpolate between linear & min phase */
@@ -278,7 +278,7 @@
work[i ] = x * cos(work[i + 1]);
work[i + 1] = x * sin(work[i + 1]);
}
- rdft(work_len, -1, work, p->bit_rev_table, p->sin_cos_table);
+ lsx_rdft(work_len, -1, work, p->bit_rev_table, p->sin_cos_table);
for (i = 0; i < work_len; ++i) work[i] *= 2. / work_len;
for (i = 1; i < work_len; ++i) if (work[i] > work[peak]) /* Find peak. */
@@ -359,7 +359,7 @@
p->bit_rev_table = calloc(dft_br_len(dft_length),sizeof(*p->bit_rev_table));
p->sin_cos_table = calloc(dft_sc_len(dft_length),sizeof(*p->sin_cos_table));
}
- rdft(dft_length, 1, f->coefs, p->bit_rev_table, p->sin_cos_table);
+ lsx_rdft(dft_length, 1, f->coefs, p->bit_rev_table, p->sin_cos_table);
}
#include "rate_filters.h"
--- a/src/sox_i.h
+++ b/src/sox_i.h
@@ -61,7 +61,19 @@
unsigned lsx_gcd(unsigned a, unsigned b);
unsigned lsx_lcm(unsigned a, unsigned b);
-double bessel_I_0(double x);
+double lsx_bessel_I_0(double x);
+extern int * lsx_fft_br;
+extern double * lsx_fft_sc;
+void lsx_safe_rdft(int len, int type, double * d);
+void lsx_safe_cdft(int len, int type, double * d);
+void lsx_power_spectrum(int n, double const * in, double * out);
+void lsx_power_spectrum_f(int n, float const * in, float * out);
+void lsx_apply_hann_f(float h[], const int num_points);
+void lsx_apply_hann(double h[], const int num_points);
+void lsx_apply_hamming(double h[], const int num_points);
+void lsx_apply_bartlett(double h[], const int num_points);
+double lsx_kaiser_beta(double att);
+void lsx_apply_kaiser(double h[], const int num_points, double beta);
#ifndef HAVE_STRCASECMP
int strcasecmp(const char *s1, const char *s2);
--- a/src/spectrogram.c
+++ b/src/spectrogram.c
@@ -40,49 +40,6 @@
#define MAX_DFT_SIZE (DFT_BASE_SIZE << MAX_DFT_SIZE_SHIFT)
#define MAX_COLS 999 /* Also max seconds */
-static void apply_hann(double h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- double x = 2 * M_PI * i / m;
- h[i] *= .5 - .5 * cos(x);
- }
-}
-
-static void apply_hamming(double h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- double x = 2 * M_PI * i / m;
- h[i] *= .53836 - .46164 * cos(x);
- }
-}
-
-static void apply_bartlett(double h[], const int num_points)
-{- int i, m = num_points - 1;
- for (i = 0; i < num_points; ++i) {- h[i] *= 2. / m * (m / 2. - fabs(i - m / 2.));
- }
-}
-
-static double kaiser_beta(double att)
-{- if (att > 100 ) return .1117 * att - 1.11;
- if (att > 50 ) return .1102 * (att - 8.7);
- if (att > 20.96) return .58417 * pow(att -20.96, .4) + .07886 * (att - 20.96);
- return 0;
-}
-
-static void apply_kaiser(double h[], const int num_points, double beta)
-{- int i, m = num_points - 1;
- for (i = 0; i <= m; ++i) {- double x = 2. * i / m - 1;
- h[i] *= bessel_I_0(beta * sqrt(1 - x * x)) / bessel_I_0(beta);
- }
-}
-
typedef enum {Window_Hann, Window_Hamming, Window_Bartlett, Window_Rectangular, Window_Kaiser} win_type_t; static enum_item const window_options[] = {ENUM_ITEM(Window_,Hann)
@@ -178,11 +135,11 @@
if (end) memset(p->window, 0, sizeof(p->window));
for (i = 0; i < n; ++i) w[i] = 1;
switch (p->win_type) {- case Window_Hann: apply_hann(w, n); break;
- case Window_Hamming: apply_hamming(w, n); break;
- case Window_Bartlett: apply_bartlett(w, n); break;
+ case Window_Hann: lsx_apply_hann(w, n); break;
+ case Window_Hamming: lsx_apply_hamming(w, n); break;
+ case Window_Bartlett: lsx_apply_bartlett(w, n); break;
case Window_Rectangular: break;
- default: apply_kaiser(w, n, kaiser_beta(p->dB_range + 20.));
+ default: lsx_apply_kaiser(w, n, lsx_kaiser_beta(p->dB_range + 20.));
}
for (i = 0; i < p->dft_size; ++i) sum += p->window[i];
for (i = 0; i < p->dft_size; ++i) p->window[i] *= 2 / sum
@@ -265,7 +222,7 @@
if ((p->end = max(p->end, p->end_min)) != p->last_end)
make_window(p, p->last_end = p->end);
for (i = 0; i < p->dft_size; ++i) p->dft_buf[i] = p->buf[i] * p->window[i];
- rdft(p->dft_size, 1, p->dft_buf, p->bit_rev_table, p->sin_cos_table);
+ lsx_rdft(p->dft_size, 1, p->dft_buf, p->bit_rev_table, p->sin_cos_table);
p->magnitudes[0] += sqr(p->dft_buf[0]);
for (i = 1; i < p->dft_size >> 1; ++i)
p->magnitudes[i] += sqr(p->dft_buf[2*i]) + sqr(p->dft_buf[2*i+1]);
--- a/src/stat.c
+++ b/src/stat.c
@@ -15,7 +15,6 @@
#include "sox_i.h"
#include <string.h>
-#include "FFT.h"
/* Private data for stat effect */
typedef struct {@@ -105,7 +104,7 @@
if (stat->fft) {stat->fft_offset = 0;
stat->re_in = lsx_malloc(sizeof(float) * stat->fft_size);
- stat->re_out = lsx_malloc(sizeof(float) * (stat->fft_size / 2));
+ stat->re_out = lsx_malloc(sizeof(float) * (stat->fft_size / 2 + 1));
}
return SOX_SUCCESS;
@@ -119,8 +118,8 @@
float ffa = rate / samples;
unsigned i;
- PowerSpectrum(samples, re_in, re_out);
- for (i = 0; i < samples / 2; i++)
+ lsx_power_spectrum_f((int)samples, re_in, re_out);
+ for (i = 0; i < samples / 2; i++) /* FIXME: should be <= samples / 2 */
fprintf(stderr, "%f %f\n", ffa * i, re_out[i]);
}