Re: [FFmpeg-devel] [PATCH 2/2] imdct15: replace the FFT with a faster PFA FFT algorithm
On 4 January 2017 at 22:45, Rostislav Pehlivanovwrote: > > > On 4 January 2017 at 14:14, Peter Barfuss wrote: > >> First off, many thanks. >> >> > +const int inv_1 = l_ptwo << ((4 - b_ptwo) & 3); >> > +const int inv_2 = 0xeeef & ((1U << b_ptwo) - 1); >> >> It would be nice to add a comment here that the expression for inv_1 >> is (2^b_ptwo)^-1 mod 15 and inv_2 is 15^-1 mod 2^b_ptwo. (A general >> PFA FFT would need to use extended Euclidean algorithm here, but >> because both cases are fixed, it simplifies to these expressions. I >> have a sketch of a proof (basically solving the relevant diophantine >> equation you get) in case anyone is nervous, though it's easy to >> verify by hand for 1 < b_ptwo < 18, which are all the cases that >> ffmpeg's power-of-two FFT currently supports). >> >> Rest of patch seems good. >> >> -Peter >> ___ >> ffmpeg-devel mailing list >> ffmpeg-devel@ffmpeg.org >> http://ffmpeg.org/mailman/listinfo/ffmpeg-devel >> > > Done > > Also I didn't like how s->exptab[20].re needed a negative sign so I > removed it > and moved the subtraction to the 5-point FFT, just so you know. > > > I'll push both patches tomorrow evening if no one else has anything to say. > Applied, Thanks ___ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
Re: [FFmpeg-devel] [PATCH 2/2] imdct15: replace the FFT with a faster PFA FFT algorithm
On 4 January 2017 at 14:14, Peter Barfusswrote: > First off, many thanks. > > > +const int inv_1 = l_ptwo << ((4 - b_ptwo) & 3); > > +const int inv_2 = 0xeeef & ((1U << b_ptwo) - 1); > > It would be nice to add a comment here that the expression for inv_1 > is (2^b_ptwo)^-1 mod 15 and inv_2 is 15^-1 mod 2^b_ptwo. (A general > PFA FFT would need to use extended Euclidean algorithm here, but > because both cases are fixed, it simplifies to these expressions. I > have a sketch of a proof (basically solving the relevant diophantine > equation you get) in case anyone is nervous, though it's easy to > verify by hand for 1 < b_ptwo < 18, which are all the cases that > ffmpeg's power-of-two FFT currently supports). > > Rest of patch seems good. > > -Peter > ___ > ffmpeg-devel mailing list > ffmpeg-devel@ffmpeg.org > http://ffmpeg.org/mailman/listinfo/ffmpeg-devel > Done Also I didn't like how s->exptab[20].re needed a negative sign so I removed it and moved the subtraction to the 5-point FFT, just so you know. I'll push both patches tomorrow evening if no one else has anything to say. ___ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
Re: [FFmpeg-devel] [PATCH 2/2] imdct15: replace the FFT with a faster PFA FFT algorithm
First off, many thanks. > +const int inv_1 = l_ptwo << ((4 - b_ptwo) & 3); > +const int inv_2 = 0xeeef & ((1U << b_ptwo) - 1); It would be nice to add a comment here that the expression for inv_1 is (2^b_ptwo)^-1 mod 15 and inv_2 is 15^-1 mod 2^b_ptwo. (A general PFA FFT would need to use extended Euclidean algorithm here, but because both cases are fixed, it simplifies to these expressions. I have a sketch of a proof (basically solving the relevant diophantine equation you get) in case anyone is nervous, though it's easy to verify by hand for 1 < b_ptwo < 18, which are all the cases that ffmpeg's power-of-two FFT currently supports). Rest of patch seems good. -Peter ___ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
Re: [FFmpeg-devel] [PATCH 2/2] imdct15: replace the FFT with a faster PFA FFT algorithm
> +/* 15-point FFT exptab */ > +for (i = 0; i < 19; i++) { > +if (i < 15) { > +double theta = 2 * M_PI * i / 15; > +s->exptab[i].re = cos(theta); > +s->exptab[i].im = sin(theta); > +} else { /* Wrap around to simplify fft[1]5 */ > +s->exptab[i] = s->exptab[i - 15]; Slightly surprised that you don't just leave these pregenerated/hardcoded as it's just 24 float values. Is this for ease of making a fixed-point version later, or some other reason? ___ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
Re: [FFmpeg-devel] [PATCH 2/2] imdct15: replace the FFT with a faster PFA FFT algorithm
On 2017-01-04 13:17, Rostislav Pehlivanov wrote: > Forgot to check the return value here, changed locally to: > > if (ff_fft_init(>ptwo_fft, N - 1, 1) < 0); > goto fail; I hope you have not changed it to that, with that semicolon at the end of the line. signature.asc Description: OpenPGP digital signature ___ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
Re: [FFmpeg-devel] [PATCH 2/2] imdct15: replace the FFT with a faster PFA FFT algorithm
On 4 January 2017 at 10:16, Rostislav Pehlivanovwrote: > > +ff_fft_init(>ptwo_fft, N - 1, 1); > > Forgot to check the return value here, changed locally to: if (ff_fft_init(>ptwo_fft, N - 1, 1) < 0); goto fail; ___ ffmpeg-devel mailing list ffmpeg-devel@ffmpeg.org http://ffmpeg.org/mailman/listinfo/ffmpeg-devel
[FFmpeg-devel] [PATCH 2/2] imdct15: replace the FFT with a faster PFA FFT algorithm
This commit replaces the current inefficient non-power-of-two FFT with a much faster FFT based on the Prime Factor Algorithm. Although it is already much faster than the old algorithm without SIMD, the new algorithm makes use of the already very throughouly SIMD'd power of two FFT, which improves performance even more across all platforms which we have SIMD support for. Most of the work was done by Peter Barfuss, who passed the code to me to implement into the iMDCT and the current codebase. The code for a 5-point and 15-point FFT was derived from the previous implementation, although it was optimized and simplified, which will make its future SIMD easier. The 15-point FFT is currently using 6% of the current overall decoder overhead. The FFT can now easily be used as a forward transform by simply not multiplying the 5-point FFT's imaginary component by -1 (which comes from the fact that changing the complex exponential's angle by -1 also changes the output by that) and by multiplying the "theta" angle of the main exptab by -1. Hence the deliberately left multiplication by -1 at the end. FATE passes, and performance reports on other platforms/CPUs are welcome. Performance comparisons: iMDCT, PFA: 101127 decicycles in speed, 32765 runs, 3 skips iMDCT, Old: 211022 decicycles in speed, 32768 runs, 0 skips Standalone FFT, 30 transforms of size 960: PFAOld FFT kiss_fftlibfftw3f 3.659695s, 15.726912s, 13.300789s, 1.18s Being only 3x slower than libfftw3f is a big achievement by itself. There appears to be something capping the performance in the iMDCT side of things, possibly during the pre-stage reindexing. However, it is certainly fast enough for now. Signed-off-by: Rostislav Pehlivanov--- libavcodec/imdct15.c | 299 ++- libavcodec/imdct15.h | 9 +- 2 files changed, 160 insertions(+), 148 deletions(-) diff --git a/libavcodec/imdct15.c b/libavcodec/imdct15.c index 7481c026cf..6d453d5cad 100644 --- a/libavcodec/imdct15.c +++ b/libavcodec/imdct15.c @@ -34,10 +34,6 @@ #include "avfft.h" #include "imdct15.h" -#include "opus.h" - -// minimal iMDCT size to make SIMD opts easier -#define CELT_MIN_IMDCT_SIZE 120 // complex c = a * b #define CMUL3(cre, cim, are, aim, bre, bim) \ @@ -48,37 +44,18 @@ do { \ #define CMUL(c, a, b) CMUL3((c).re, (c).im, (a).re, (a).im, (b).re, (b).im) -// complex c = a * b -// d = a * conjugate(b) -#define CMUL2(c, d, a, b)\ -do { \ -float are = (a).re; \ -float aim = (a).im; \ -float bre = (b).re; \ -float bim = (b).im; \ -float rr = are * bre; \ -float ri = are * bim; \ -float ir = aim * bre; \ -float ii = aim * bim; \ -(c).re = rr - ii; \ -(c).im = ri + ir; \ -(d).re = rr + ii; \ -(d).im = -ri + ir; \ -} while (0) - av_cold void ff_imdct15_uninit(IMDCT15Context **ps) { IMDCT15Context *s = *ps; -int i; if (!s) return; -for (i = 0; i < FF_ARRAY_ELEMS(s->exptab); i++) -av_freep(>exptab[i]); +ff_fft_end(>ptwo_fft); +av_freep(>pfa_prereindex); +av_freep(>pfa_postreindex); av_freep(>twiddle_exptab); - av_freep(>tmp); av_freep(ps); @@ -87,14 +64,46 @@ av_cold void ff_imdct15_uninit(IMDCT15Context **ps) static void imdct15_half(IMDCT15Context *s, float *dst, const float *src, ptrdiff_t stride, float scale); +static inline int init_pfa_reindex_tabs(IMDCT15Context *s) +{ +int i, j; +const int b_ptwo = s->ptwo_fft.nbits; /* Bits for the power of two FFTs */ +const int l_ptwo = 1 << b_ptwo; /* Total length for the power of two FFTs */ +const int inv_1 = l_ptwo << ((4 - b_ptwo) & 3); +const int inv_2 = 0xeeef & ((1U << b_ptwo) - 1); + +s->pfa_prereindex = av_malloc(15 * l_ptwo * sizeof(*s->pfa_prereindex)); +if (!s->pfa_prereindex) +return 1; + +s->pfa_postreindex = av_malloc(15 * l_ptwo * sizeof(*s->pfa_postreindex)); +if (!s->pfa_postreindex) +return 1; + +/* Pre/Post-reindex */ +for (i = 0; i < l_ptwo; i++) { +for (j = 0; j < 15; j++) { +const int q_pre = ((l_ptwo * j)/15 + i) >> b_ptwo; +const int q_post = (((j*inv_1)/15) + (i*inv_2)) >> b_ptwo; +const int k_pre = 15*i + ((j - q_pre*15) << b_ptwo); +const int k_post = i*inv_2*15 + j*inv_1 - 15*q_post*l_ptwo; +s->pfa_prereindex [i*15 + j] = k_pre; +