lavu/tx: add support for double precision FFT and MDCT

Simply moves and templates the actual transforms to support an additional data type. Unlike the float version, which is equal or better than libfftw3f, double precision output is bit identical with libfftw3.

Lynne committed 7 years ago

Commit 42e23191 parent f60b1219 files · +825 −677Browse files

Unified · Split

9 files

doc/APIchanges +3 −0
libavutil/Makefile +2 −0
libavutil/tx.c +17 −676
libavutil/tx.h +12 −0
libavutil/tx_double.c +21 −0
libavutil/tx_float.c +21 −0
libavutil/tx_priv.h +105 −0
libavutil/tx_template.c +643 −0
libavutil/version.h +1 −1

doc/APIchanges+3−0

@@ -15,6 +15,9 @@ libavutil: 2017-10-21 · 14 unmodified lines
15	15
16	16	API changes, most recent first:
17	17
	18	+2019-07-27 - xxxxxxxxxx - lavu 56.33.100 - tx.h
	19	+ Add AV_TX_DOUBLE_FFT and AV_TX_DOUBLE_MDCT
	20	+
18	21	-------- 8< --------- FFmpeg 4.2 was cut here -------- 8< ---------
19	22
20	23	2019-06-21 - a30e44098a - lavu 56.30.100 - frame.h

libavutil/Makefile+2−0

@@ -161,6 +161,8 @@ OBJS = adler32.o \ · 160 unmodified lines
161	161	xtea.o \
162	162	tea.o \
163	163	tx.o \
	164	+ tx_float.o \
	165	+ tx_double.o
164	166
165	167	OBJS-$(CONFIG_CUDA) += hwcontext_cuda.o
166	168	OBJS-$(CONFIG_D3D11VA) += hwcontext_d3d11va.o

libavutil/tx.c+17−676

@@ -1,10 +1,4 @@
1	1	/*
2		− * Copyright (c) 2019 Lynne <dev@lynne.ee>
3		− * Power of two FFT:
4		− * Copyright (c) 2008 Loren Merritt
5		− * Copyright (c) 2002 Fabrice Bellard
6		− * Partly based on libdjbfft by D. J. Bernstein
7		− *
8	2	* This file is part of FFmpeg.
9	3	*
10	4	* FFmpeg is free software; you can redistribute it and/or
@@ -22,576 +16,10 @@ · 11 unmodified lines
22	16	* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23	17	*/
24	18
25		−#include <stddef.h>
26		−#include "tx.h"
27		−#include "thread.h"
28		−#include "mem.h"
29		−#include "avassert.h"
30		−
31		−typedef float FFTSample;
32		−typedef AVComplexFloat FFTComplex;
33		−
34		−struct AVTXContext {
35		− int n; /* Nptwo part */
36		− int m; /* Ptwo part */
37		−
38		− FFTComplex exptab; / MDCT exptab */
39		− FFTComplex tmp; / Temporary buffer needed for all compound transforms */
40		− int pfatab; / Input/Output mapping for compound transforms */
41		− int revtab; / Input mapping for power of two transforms */
42		−};
43		−
44		−#define FFT_NAME(x) x
45		−
46		−#define COSTABLE(size) \
47		− static DECLARE_ALIGNED(32, FFTSample, FFT_NAME(ff_cos_##size))[size/2]
48		−
49		−static FFTSample * const FFT_NAME(ff_cos_tabs)[18];
50		−
51		−COSTABLE(16);
52		−COSTABLE(32);
53		−COSTABLE(64);
54		−COSTABLE(128);
55		−COSTABLE(256);
56		−COSTABLE(512);
57		−COSTABLE(1024);
58		−COSTABLE(2048);
59		−COSTABLE(4096);
60		−COSTABLE(8192);
61		−COSTABLE(16384);
62		−COSTABLE(32768);
63		−COSTABLE(65536);
64		−COSTABLE(131072);
65		−
66		−static av_cold void init_ff_cos_tabs(int index)
67		−{
68		− int m = 1 << index;
69		− double freq = 2*M_PI/m;
70		− FFTSample *tab = FFT_NAME(ff_cos_tabs)[index];
71		− for(int i = 0; i <= m/4; i++)
72		− tab[i] = cos(i*freq);
73		− for(int i = 1; i < m/4; i++)
74		− tab[m/2 - i] = tab[i];
75		−}
76		−
77		−typedef struct CosTabsInitOnce {
78		− void (*func)(void);
79		− AVOnce control;
80		−} CosTabsInitOnce;
81		−
82		−#define INIT_FF_COS_TABS_FUNC(index, size) \
83		−static av_cold void init_ff_cos_tabs_ ## size (void) \
84		−{ \
85		− init_ff_cos_tabs(index); \
86		−}
87		−
88		−INIT_FF_COS_TABS_FUNC(4, 16)
89		−INIT_FF_COS_TABS_FUNC(5, 32)
90		−INIT_FF_COS_TABS_FUNC(6, 64)
91		−INIT_FF_COS_TABS_FUNC(7, 128)
92		−INIT_FF_COS_TABS_FUNC(8, 256)
93		−INIT_FF_COS_TABS_FUNC(9, 512)
94		−INIT_FF_COS_TABS_FUNC(10, 1024)
95		−INIT_FF_COS_TABS_FUNC(11, 2048)
96		−INIT_FF_COS_TABS_FUNC(12, 4096)
97		−INIT_FF_COS_TABS_FUNC(13, 8192)
98		−INIT_FF_COS_TABS_FUNC(14, 16384)
99		−INIT_FF_COS_TABS_FUNC(15, 32768)
100		−INIT_FF_COS_TABS_FUNC(16, 65536)
101		−INIT_FF_COS_TABS_FUNC(17, 131072)
102		−
103		−static CosTabsInitOnce cos_tabs_init_once[] = {
104		− { NULL },
105		− { NULL },
106		− { NULL },
107		− { NULL },
108		− { init_ff_cos_tabs_16, AV_ONCE_INIT },
109		− { init_ff_cos_tabs_32, AV_ONCE_INIT },
110		− { init_ff_cos_tabs_64, AV_ONCE_INIT },
111		− { init_ff_cos_tabs_128, AV_ONCE_INIT },
112		− { init_ff_cos_tabs_256, AV_ONCE_INIT },
113		− { init_ff_cos_tabs_512, AV_ONCE_INIT },
114		− { init_ff_cos_tabs_1024, AV_ONCE_INIT },
115		− { init_ff_cos_tabs_2048, AV_ONCE_INIT },
116		− { init_ff_cos_tabs_4096, AV_ONCE_INIT },
117		− { init_ff_cos_tabs_8192, AV_ONCE_INIT },
118		− { init_ff_cos_tabs_16384, AV_ONCE_INIT },
119		− { init_ff_cos_tabs_32768, AV_ONCE_INIT },
120		− { init_ff_cos_tabs_65536, AV_ONCE_INIT },
121		− { init_ff_cos_tabs_131072, AV_ONCE_INIT },
122		−};
123		−
124		−static FFTSample * const FFT_NAME(ff_cos_tabs)[] = {
125		− NULL, NULL, NULL, NULL,
126		− FFT_NAME(ff_cos_16),
127		− FFT_NAME(ff_cos_32),
128		− FFT_NAME(ff_cos_64),
129		− FFT_NAME(ff_cos_128),
130		− FFT_NAME(ff_cos_256),
131		− FFT_NAME(ff_cos_512),
132		− FFT_NAME(ff_cos_1024),
133		− FFT_NAME(ff_cos_2048),
134		− FFT_NAME(ff_cos_4096),
135		− FFT_NAME(ff_cos_8192),
136		− FFT_NAME(ff_cos_16384),
137		− FFT_NAME(ff_cos_32768),
138		− FFT_NAME(ff_cos_65536),
139		− FFT_NAME(ff_cos_131072),
140		−};
141		−
142		−static av_cold void ff_init_ff_cos_tabs(int index)
143		−{
144		− ff_thread_once(&cos_tabs_init_once[index].control,
145		− cos_tabs_init_once[index].func);
146		−}
147		−
148		−static AVOnce tabs_53_once = AV_ONCE_INIT;
149		−static DECLARE_ALIGNED(32, FFTComplex, FFT_NAME(ff_53_tabs))[4];
150		−
151		−static av_cold void ff_init_53_tabs(void)
152		−{
153		− ff_53_tabs[0] = (FFTComplex){ cos(2 * M_PI / 12), cos(2 * M_PI / 12) };
154		− ff_53_tabs[1] = (FFTComplex){ 0.5, 0.5 };
155		− ff_53_tabs[2] = (FFTComplex){ cos(2 * M_PI / 5), sin(2 * M_PI / 5) };
156		− ff_53_tabs[3] = (FFTComplex){ cos(2 * M_PI / 10), sin(2 * M_PI / 10) };
157		−}
158		−
159		−#define BF(x, y, a, b) do { \
160		− x = (a) - (b); \
161		− y = (a) + (b); \
162		− } while (0)
163		−
164		−#define CMUL(dre, dim, are, aim, bre, bim) do { \
165		− (dre) = (are) * (bre) - (aim) * (bim); \
166		− (dim) = (are) * (bim) + (aim) * (bre); \
167		− } while (0)
168		−
169		−#define CMUL3(c, a, b) CMUL((c).re, (c).im, (a).re, (a).im, (b).re, (b).im)
170		−
171		−static av_always_inline void fft3(FFTComplex out, FFTComplex in,
172		− ptrdiff_t stride)
173		−{
174		− FFTComplex tmp[2];
175		−
176		− tmp[0].re = in[1].im - in[2].im;
177		− tmp[0].im = in[1].re - in[2].re;
178		− tmp[1].re = in[1].re + in[2].re;
179		− tmp[1].im = in[1].im + in[2].im;
180		−
181		− out[0*stride].re = in[0].re + tmp[1].re;
182		− out[0*stride].im = in[0].im + tmp[1].im;
183		−
184		− tmp[0].re *= ff_53_tabs[0].re;
185		− tmp[0].im *= ff_53_tabs[0].im;
186		− tmp[1].re *= ff_53_tabs[1].re;
187		− tmp[1].im *= ff_53_tabs[1].re;
188		−
189		− out[1*stride].re = in[0].re - tmp[1].re + tmp[0].re;
190		− out[1*stride].im = in[0].im - tmp[1].im - tmp[0].im;
191		− out[2*stride].re = in[0].re - tmp[1].re - tmp[0].re;
192		− out[2*stride].im = in[0].im - tmp[1].im + tmp[0].im;
193		−}
194		−
195		−#define DECL_FFT5(NAME, D0, D1, D2, D3, D4) \
196		−static av_always_inline void NAME(FFTComplex out, FFTComplex in, \
197		− ptrdiff_t stride) \
198		−{ \
199		− FFTComplex z0[4], t[6]; \
200		− \
201		− t[0].re = in[1].re + in[4].re; \
202		− t[0].im = in[1].im + in[4].im; \
203		− t[1].im = in[1].re - in[4].re; \
204		− t[1].re = in[1].im - in[4].im; \
205		− t[2].re = in[2].re + in[3].re; \
206		− t[2].im = in[2].im + in[3].im; \
207		− t[3].im = in[2].re - in[3].re; \
208		− t[3].re = in[2].im - in[3].im; \
209		− \
210		− out[D0*stride].re = in[0].re + in[1].re + in[2].re + \
211		− in[3].re + in[4].re; \
212		− out[D0*stride].im = in[0].im + in[1].im + in[2].im + \
213		− in[3].im + in[4].im; \
214		− \
215		− t[4].re = ff_53_tabs[2].re * t[2].re - ff_53_tabs[3].re * t[0].re; \
216		− t[4].im = ff_53_tabs[2].re * t[2].im - ff_53_tabs[3].re * t[0].im; \
217		− t[0].re = ff_53_tabs[2].re * t[0].re - ff_53_tabs[3].re * t[2].re; \
218		− t[0].im = ff_53_tabs[2].re * t[0].im - ff_53_tabs[3].re * t[2].im; \
219		− t[5].re = ff_53_tabs[2].im * t[3].re - ff_53_tabs[3].im * t[1].re; \
220		− t[5].im = ff_53_tabs[2].im * t[3].im - ff_53_tabs[3].im * t[1].im; \
221		− t[1].re = ff_53_tabs[2].im * t[1].re + ff_53_tabs[3].im * t[3].re; \
222		− t[1].im = ff_53_tabs[2].im * t[1].im + ff_53_tabs[3].im * t[3].im; \
223		− \
224		− z0[0].re = t[0].re - t[1].re; \
225		− z0[0].im = t[0].im - t[1].im; \
226		− z0[1].re = t[4].re + t[5].re; \
227		− z0[1].im = t[4].im + t[5].im; \
228		− \
229		− z0[2].re = t[4].re - t[5].re; \
230		− z0[2].im = t[4].im - t[5].im; \
231		− z0[3].re = t[0].re + t[1].re; \
232		− z0[3].im = t[0].im + t[1].im; \
233		− \
234		− out[D1*stride].re = in[0].re + z0[3].re; \
235		− out[D1*stride].im = in[0].im + z0[0].im; \
236		− out[D2*stride].re = in[0].re + z0[2].re; \
237		− out[D2*stride].im = in[0].im + z0[1].im; \
238		− out[D3*stride].re = in[0].re + z0[1].re; \
239		− out[D3*stride].im = in[0].im + z0[2].im; \
240		− out[D4*stride].re = in[0].re + z0[0].re; \
241		− out[D4*stride].im = in[0].im + z0[3].im; \
242		−}
243		−
244		−DECL_FFT5(fft5, 0, 1, 2, 3, 4)
245		−DECL_FFT5(fft5_m1, 0, 6, 12, 3, 9)
246		−DECL_FFT5(fft5_m2, 10, 1, 7, 13, 4)
247		−DECL_FFT5(fft5_m3, 5, 11, 2, 8, 14)
248		−
249		−static av_always_inline void fft15(FFTComplex out, FFTComplex in,
250		− ptrdiff_t stride)
251		−{
252		− FFTComplex tmp[15];
253		−
254		− for (int i = 0; i < 5; i++)
255		− fft3(tmp + i, in + i*3, 5);
256		−
257		− fft5_m1(out, tmp + 0, stride);
258		− fft5_m2(out, tmp + 5, stride);
259		− fft5_m3(out, tmp + 10, stride);
260		−}
261		−
262		−#define BUTTERFLIES(a0,a1,a2,a3) {\
263		− BF(t3, t5, t5, t1);\
264		− BF(a2.re, a0.re, a0.re, t5);\
265		− BF(a3.im, a1.im, a1.im, t3);\
266		− BF(t4, t6, t2, t6);\
267		− BF(a3.re, a1.re, a1.re, t4);\
268		− BF(a2.im, a0.im, a0.im, t6);\
269		−}
270		−
271		−// force loading all the inputs before storing any.
272		−// this is slightly slower for small data, but avoids store->load aliasing
273		−// for addresses separated by large powers of 2.
274		−#define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
275		− FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
276		− BF(t3, t5, t5, t1);\
277		− BF(a2.re, a0.re, r0, t5);\
278		− BF(a3.im, a1.im, i1, t3);\
279		− BF(t4, t6, t2, t6);\
280		− BF(a3.re, a1.re, r1, t4);\
281		− BF(a2.im, a0.im, i0, t6);\
282		−}
283		−
284		−#define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
285		− CMUL(t1, t2, a2.re, a2.im, wre, -wim);\
286		− CMUL(t5, t6, a3.re, a3.im, wre, wim);\
287		− BUTTERFLIES(a0,a1,a2,a3)\
288		−}
289		−
290		−#define TRANSFORM_ZERO(a0,a1,a2,a3) {\
291		− t1 = a2.re;\
292		− t2 = a2.im;\
293		− t5 = a3.re;\
294		− t6 = a3.im;\
295		− BUTTERFLIES(a0,a1,a2,a3)\
296		−}
297		−
298		−/* z[0...8n-1], w[1...2n-1] */
299		−#define PASS(name)\
300		−static void name(FFTComplex z, const FFTSample wre, unsigned int n)\
301		−{\
302		− FFTSample t1, t2, t3, t4, t5, t6;\
303		− int o1 = 2*n;\
304		− int o2 = 4*n;\
305		− int o3 = 6*n;\
306		− const FFTSample *wim = wre+o1;\
307		− n--;\
308		−\
309		− TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
310		− TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
311		− do {\
312		− z += 2;\
313		− wre += 2;\
314		− wim -= 2;\
315		− TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
316		− TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
317		− } while(--n);\
318		−}
319		−
320		−PASS(pass)
321		−#undef BUTTERFLIES
322		−#define BUTTERFLIES BUTTERFLIES_BIG
323		−PASS(pass_big)
324		−
325		−#define DECL_FFT(n,n2,n4)\
326		−static void fft##n(FFTComplex *z)\
327		−{\
328		− fft##n2(z);\
329		− fft##n4(z+n4*2);\
330		− fft##n4(z+n4*3);\
331		− pass(z,FFT_NAME(ff_cos_##n),n4/2);\
332		−}
333		−
334		−static void fft4(FFTComplex *z)
335		−{
336		− FFTSample t1, t2, t3, t4, t5, t6, t7, t8;
337		−
338		− BF(t3, t1, z[0].re, z[1].re);
339		− BF(t8, t6, z[3].re, z[2].re);
340		− BF(z[2].re, z[0].re, t1, t6);
341		− BF(t4, t2, z[0].im, z[1].im);
342		− BF(t7, t5, z[2].im, z[3].im);
343		− BF(z[3].im, z[1].im, t4, t8);
344		− BF(z[3].re, z[1].re, t3, t7);
345		− BF(z[2].im, z[0].im, t2, t5);
346		−}
347		−
348		−static void fft8(FFTComplex *z)
349		−{
350		− FFTSample t1, t2, t3, t4, t5, t6;
351		−
352		− fft4(z);
353		−
354		− BF(t1, z[5].re, z[4].re, -z[5].re);
355		− BF(t2, z[5].im, z[4].im, -z[5].im);
356		− BF(t5, z[7].re, z[6].re, -z[7].re);
357		− BF(t6, z[7].im, z[6].im, -z[7].im);
358		−
359		− BUTTERFLIES(z[0],z[2],z[4],z[6]);
360		− TRANSFORM(z[1],z[3],z[5],z[7],M_SQRT1_2,M_SQRT1_2);
361		−}
362		−
363		−static void fft16(FFTComplex *z)
364		−{
365		− FFTSample t1, t2, t3, t4, t5, t6;
366		− FFTSample cos_16_1 = FFT_NAME(ff_cos_16)[1];
367		− FFTSample cos_16_3 = FFT_NAME(ff_cos_16)[3];
368		−
369		− fft8(z);
370		− fft4(z+8);
371		− fft4(z+12);
372		−
373		− TRANSFORM_ZERO(z[0],z[4],z[8],z[12]);
374		− TRANSFORM(z[2],z[6],z[10],z[14],M_SQRT1_2,M_SQRT1_2);
375		− TRANSFORM(z[1],z[5],z[9],z[13],cos_16_1,cos_16_3);
376		− TRANSFORM(z[3],z[7],z[11],z[15],cos_16_3,cos_16_1);
377		−}
378		−
379		−DECL_FFT(32,16,8)
380		−DECL_FFT(64,32,16)
381		−DECL_FFT(128,64,32)
382		−DECL_FFT(256,128,64)
383		−DECL_FFT(512,256,128)
384		−#define pass pass_big
385		−DECL_FFT(1024,512,256)
386		−DECL_FFT(2048,1024,512)
387		−DECL_FFT(4096,2048,1024)
388		−DECL_FFT(8192,4096,2048)
389		−DECL_FFT(16384,8192,4096)
390		−DECL_FFT(32768,16384,8192)
391		−DECL_FFT(65536,32768,16384)
392		−DECL_FFT(131072,65536,32768)
393		−
394		−static void (* const fft_dispatch[])(FFTComplex*) = {
395		− fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024,
396		− fft2048, fft4096, fft8192, fft16384, fft32768, fft65536, fft131072
397		−};
398		−
399		−#define DECL_COMP_FFT(N) \
400		−static void compound_fft_##N##xM(AVTXContext s, void _out, \
401		− void *_in, ptrdiff_t stride) \
402		−{ \
403		− const int m = s->m, in_map = s->pfatab, out_map = in_map + N*m; \
404		− FFTComplex *in = _in; \
405		− FFTComplex *out = _out; \
406		− FFTComplex fft##N##in[N]; \
407		− void (fftp)(FFTComplex z) = fft_dispatch[av_log2(m) - 2]; \
408		− \
409		− for (int i = 0; i < m; i++) { \
410		− for (int j = 0; j < N; j++) \
411		− fft##N##in[j] = in[in_map[i*N + j]]; \
412		− fft##N(s->tmp + s->revtab[i], fft##N##in, m); \
413		− } \
414		− \
415		− for (int i = 0; i < N; i++) \
416		− fftp(s->tmp + m*i); \
417		− \
418		− for (int i = 0; i < N*m; i++) \
419		− out[i] = s->tmp[out_map[i]]; \
420		−}
421		−
422		−DECL_COMP_FFT(3)
423		−DECL_COMP_FFT(5)
424		−DECL_COMP_FFT(15)
425		−
426		−static void monolithic_fft(AVTXContext s, void _out, void *_in,
427		− ptrdiff_t stride)
428		−{
429		− FFTComplex *in = _in;
430		− FFTComplex *out = _out;
431		− int m = s->m, mb = av_log2(m) - 2;
432		− for (int i = 0; i < m; i++)
433		− out[s->revtab[i]] = in[i];
434		− fft_dispatch[mb](out);
435		−}
436		−
437		−#define DECL_COMP_IMDCT(N) \
438		−static void compound_imdct_##N##xM(AVTXContext s, void _dst, void *_src, \
439		− ptrdiff_t stride) \
440		−{ \
441		− FFTComplex fft##N##in[N]; \
442		− FFTComplex z = _dst, exp = s->exptab; \
443		− const int m = s->m, len8 = N*m >> 1; \
444		− const int in_map = s->pfatab, out_map = in_map + N*m; \
445		− const float src = _src, in1, *in2; \
446		− void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2]; \
447		− \
448		− stride /= sizeof(src); / To convert it from bytes */ \
449		− in1 = src; \
450		− in2 = src + ((Nm2) - 1) * stride; \
451		− \
452		− for (int i = 0; i < m; i++) { \
453		− for (int j = 0; j < N; j++) { \
454		− const int k = in_map[i*N + j]; \
455		− FFTComplex tmp = { in2[-kstride], in1[kstride] }; \
456		− CMUL3(fft##N##in[j], tmp, exp[k >> 1]); \
457		− } \
458		− fft##N(s->tmp + s->revtab[i], fft##N##in, m); \
459		− } \
460		− \
461		− for (int i = 0; i < N; i++) \
462		− fftp(s->tmp + m*i); \
463		− \
464		− for (int i = 0; i < len8; i++) { \
465		− const int i0 = len8 + i, i1 = len8 - i - 1; \
466		− const int s0 = out_map[i0], s1 = out_map[i1]; \
467		− FFTComplex src1 = { s->tmp[s1].im, s->tmp[s1].re }; \
468		− FFTComplex src0 = { s->tmp[s0].im, s->tmp[s0].re }; \
469		− \
470		− CMUL(z[i1].re, z[i0].im, src1.re, src1.im, exp[i1].im, exp[i1].re); \
471		− CMUL(z[i0].re, z[i1].im, src0.re, src0.im, exp[i0].im, exp[i0].re); \
472		− } \
473		−}
474		−
475		−DECL_COMP_IMDCT(3)
476		−DECL_COMP_IMDCT(5)
477		−DECL_COMP_IMDCT(15)
478		−
479		−#define DECL_COMP_MDCT(N) \
480		−static void compound_mdct_##N##xM(AVTXContext s, void _dst, void *_src, \
481		− ptrdiff_t stride) \
482		−{ \
483		− float src = _src, dst = _dst; \
484		− FFTComplex *exp = s->exptab, tmp, fft##N##in[N]; \
485		− const int m = s->m, len4 = Nm, len3 = len4 3, len8 = len4 >> 1; \
486		− const int in_map = s->pfatab, out_map = in_map + N*m; \
487		− void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2]; \
488		− \
489		− stride /= sizeof(*dst); \
490		− \
491		− for (int i = 0; i < m; i++) { /* Folding and pre-reindexing */ \
492		− for (int j = 0; j < N; j++) { \
493		− const int k = in_map[i*N + j]; \
494		− if (k < len4) { \
495		− tmp.re = -src[ len4 + k] + src[1*len4 - 1 - k]; \
496		− tmp.im = -src[ len3 + k] - src[1*len3 - 1 - k]; \
497		− } else { \
498		− tmp.re = -src[ len4 + k] - src[5*len4 - 1 - k]; \
499		− tmp.im = src[-len4 + k] - src[1*len3 - 1 - k]; \
500		− } \
501		− CMUL(fft##N##in[j].im, fft##N##in[j].re, tmp.re, tmp.im, \
502		− exp[k >> 1].re, exp[k >> 1].im); \
503		− } \
504		− fft##N(s->tmp + s->revtab[i], fft##N##in, m); \
505		− } \
506		− \
507		− for (int i = 0; i < N; i++) \
508		− fftp(s->tmp + m*i); \
509		− \
510		− for (int i = 0; i < len8; i++) { \
511		− const int i0 = len8 + i, i1 = len8 - i - 1; \
512		− const int s0 = out_map[i0], s1 = out_map[i1]; \
513		− FFTComplex src1 = { s->tmp[s1].re, s->tmp[s1].im }; \
514		− FFTComplex src0 = { s->tmp[s0].re, s->tmp[s0].im }; \
515		− \
516		− CMUL(dst[2i1stride + stride], dst[2i0stride], src0.re, src0.im, \
517		− exp[i0].im, exp[i0].re); \
518		− CMUL(dst[2i0stride + stride], dst[2i1stride], src1.re, src1.im, \
519		− exp[i1].im, exp[i1].re); \
520		− } \
521		−}
522		−
523		−DECL_COMP_MDCT(3)
524		−DECL_COMP_MDCT(5)
525		−DECL_COMP_MDCT(15)
526		−
527		−static void monolithic_imdct(AVTXContext s, void _dst, void *_src,
528		− ptrdiff_t stride)
529		−{
530		− FFTComplex z = _dst, exp = s->exptab;
531		− const int m = s->m, len8 = m >> 1;
532		− const float src = _src, in1, *in2;
533		− void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2];
534		−
535		− stride /= sizeof(*src);
536		− in1 = src;
537		− in2 = src + ((m2) - 1) stride;
538		−
539		− for (int i = 0; i < m; i++) {
540		− FFTComplex tmp = { in2[-2istride], in1[2istride] };
541		− CMUL3(z[s->revtab[i]], tmp, exp[i]);
542		− }
543		−
544		− fftp(z);
545		−
546		− for (int i = 0; i < len8; i++) {
547		− const int i0 = len8 + i, i1 = len8 - i - 1;
548		− FFTComplex src1 = { z[i1].im, z[i1].re };
549		− FFTComplex src0 = { z[i0].im, z[i0].re };
550		−
551		− CMUL(z[i1].re, z[i0].im, src1.re, src1.im, exp[i1].im, exp[i1].re);
552		− CMUL(z[i0].re, z[i1].im, src0.re, src0.im, exp[i0].im, exp[i0].re);
553		− }
554		−}
555		−
556		−static void monolithic_mdct(AVTXContext s, void _dst, void *_src,
557		− ptrdiff_t stride)
558		−{
559		− float src = _src, dst = _dst;
560		− FFTComplex exp = s->exptab, tmp, z = _dst;
561		− const int m = s->m, len4 = m, len3 = len4 * 3, len8 = len4 >> 1;
562		− void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2];
563		−
564		− stride /= sizeof(*dst);
565		−
566		− for (int i = 0; i < m; i++) { /* Folding and pre-reindexing */
567		− const int k = 2*i;
568		− if (k < len4) {
569		− tmp.re = -src[ len4 + k] + src[1*len4 - 1 - k];
570		− tmp.im = -src[ len3 + k] - src[1*len3 - 1 - k];
571		− } else {
572		− tmp.re = -src[ len4 + k] - src[5*len4 - 1 - k];
573		− tmp.im = src[-len4 + k] - src[1*len3 - 1 - k];
574		− }
575		− CMUL(z[s->revtab[i]].im, z[s->revtab[i]].re, tmp.re, tmp.im,
576		− exp[i].re, exp[i].im);
577		− }
578		−
579		− fftp(z);
580		−
581		− for (int i = 0; i < len8; i++) {
582		− const int i0 = len8 + i, i1 = len8 - i - 1;
583		− FFTComplex src1 = { z[i1].re, z[i1].im };
584		− FFTComplex src0 = { z[i0].re, z[i0].im };
585		−
586		− CMUL(dst[2i1stride + stride], dst[2i0stride], src0.re, src0.im,
587		− exp[i0].im, exp[i0].re);
588		− CMUL(dst[2i0stride + stride], dst[2i1stride], src1.re, src1.im,
589		− exp[i1].im, exp[i1].re);
590		− }
591		−}
	19	+#include "tx_priv.h"
592	20
593	21	/* Calculates the modular multiplicative inverse, not fast, replace */
594		−static int mulinv(int n, int m)
	22	+static av_always_inline int mulinv(int n, int m)
595	23	{
596	24	n = n % m;
597	25	for (int x = 1; x < m; x++)
@@ -601,14 +29,17 @@ static int mulinv(int n, int m) · 3 unmodified lines
601	29	}
602	30
603	31	/* Guaranteed to work for any n, m where gcd(n, m) == 1 */
604		−static int gen_compound_mapping(AVTXContext *s, int n, int m, int inv,
605		− enum AVTXType type)
	32	+int ff_tx_gen_compound_mapping(AVTXContext *s)
606	33	{
607	34	int in_map, out_map;
	35	+ const int n = s->n;
	36	+ const int m = s->m;
	37	+ const int inv = s->inv;
	38	+ const int type = s->type;
608	39	const int len = n*m;
609	40	const int m_inv = mulinv(m, n);
610	41	const int n_inv = mulinv(n, m);
611		− const int mdct = type == AV_TX_FLOAT_MDCT;
	42	+ const int mdct = type == AV_TX_FLOAT_MDCT \|\| type == AV_TX_DOUBLE_MDCT;
612	43
613	44	if (!(s->pfatab = av_malloc(2lensizeof(*s->pfatab))))
614	45	return AVERROR(ENOMEM);
@@ -619,7 +50,7 @@ static int gen_compound_mapping(AVTXContext *s, int n, int m, int inv, · 4 unmodified lines
619	50	/* Ruritanian map for input, CRT map for output, can be swapped */
620	51	for (int j = 0; j < m; j++) {
621	52	for (int i = 0; i < n; i++) {
622		− /* Shifted by 1 to simplify forward MDCTs */
	53	+ /* Shifted by 1 to simplify MDCTs */
623	54	in_map[jn + i] = ((im + j*n) % len) << mdct;
624	55	out_map[(imm_inv + jnn_inv) % len] = i*m + j;
625	56	}
@@ -649,23 +80,10 @@ static int gen_compound_mapping(AVTXContext *s, int n, int m, int inv, · 23 unmodified lines
649	80	return 0;
650	81	}
651	82
652		−static int split_radix_permutation(int i, int n, int inverse)
	83	+int ff_tx_gen_ptwo_revtab(AVTXContext *s)
653	84	{
654		− int m;
655		− if (n <= 2)
656		− return i & 1;
657		− m = n >> 1;
658		− if (!(i & m))
659		− return split_radix_permutation(i, m, inverse)*2;
660		− m >>= 1;
661		− if (inverse == !(i & m))
662		− return split_radix_permutation(i, m, inverse)*4 + 1;
663		− else
664		− return split_radix_permutation(i, m, inverse)*4 - 1;
665		−}
	85	+ const int m = s->m, inv = s->inv;
666	86
667		−static int get_ptwo_revtab(AVTXContext *s, int m, int inv)
668		−{
669	87	if (!(s->revtab = av_malloc(msizeof(s->revtab))))
670	88	return AVERROR(ENOMEM);
671	89
@@ -678,23 +96,6 @@ static int get_ptwo_revtab(AVTXContext *s, int m, int inv) · 6 unmodified lines
678	96	return 0;
679	97	}
680	98
681		−static int gen_mdct_exptab(AVTXContext *s, int len4, double scale)
682		−{
683		− const double theta = (scale < 0 ? len4 : 0) + 1.0/8.0;
684		−
685		− if (!(s->exptab = av_malloc_array(len4, sizeof(*s->exptab))))
686		− return AVERROR(ENOMEM);
687		−
688		− scale = sqrt(fabs(scale));
689		− for (int i = 0; i < len4; i++) {
690		− const double alpha = M_PI_2 * (i + theta) / len4;
691		− s->exptab[i].re = cos(alpha) * scale;
692		− s->exptab[i].im = sin(alpha) * scale;
693		− }
694		−
695		− return 0;
696		−}
697		−
698	99	av_cold void av_tx_uninit(AVTXContext **ctx)
699	100	{
700	101	if (!(*ctx))
@@ -708,71 +109,6 @@ av_cold void av_tx_uninit(AVTXContext **ctx) · 7 unmodified lines
708	109	av_freep(ctx);
709	110	}
710	111
711		−static int init_mdct_fft(AVTXContext s, av_tx_fn tx, enum AVTXType type,
712		− int inv, int len, const void *scale, uint64_t flags)
713		−{
714		− int err, n = 1, m = 1, max_ptwo = 1 << (FF_ARRAY_ELEMS(fft_dispatch) + 1);
715		−
716		− if (type == AV_TX_FLOAT_MDCT)
717		− len >>= 1;
718		−
719		−#define CHECK_FACTOR(DST, FACTOR, SRC) \
720		− if (DST == 1 && !(SRC % FACTOR)) { \
721		− DST = FACTOR; \
722		− SRC /= FACTOR; \
723		− }
724		− CHECK_FACTOR(n, 15, len)
725		− CHECK_FACTOR(n, 5, len)
726		− CHECK_FACTOR(n, 3, len)
727		−#undef CHECK_NPTWO_FACTOR
728		−
729		− /* len must be a power of two now */
730		− if (!(len & (len - 1)) && len >= 4 && len <= max_ptwo) {
731		− m = len;
732		− len = 1;
733		− }
734		−
735		− /* Filter out direct 3, 5 and 15 transforms, too niche */
736		− if (len > 1 \|\| m == 1) {
737		− av_log(NULL, AV_LOG_ERROR, "Unsupported transform size: n = %i, "
738		− "m = %i, residual = %i!\n", n, m, len);
739		− return AVERROR(EINVAL);
740		− } else if (n > 1 && m > 1) { /* 2D transform case */
741		− if ((err = gen_compound_mapping(s, n, m, inv, type)))
742		− return err;
743		− if (!(s->tmp = av_malloc(nmsizeof(*s->tmp))))
744		− return AVERROR(ENOMEM);
745		− *tx = n == 3 ? compound_fft_3xM :
746		− n == 5 ? compound_fft_5xM :
747		− compound_fft_15xM;
748		− if (type == AV_TX_FLOAT_MDCT)
749		− *tx = n == 3 ? inv ? compound_imdct_3xM : compound_mdct_3xM :
750		− n == 5 ? inv ? compound_imdct_5xM : compound_mdct_5xM :
751		− inv ? compound_imdct_15xM : compound_mdct_15xM;
752		− } else { /* Direct transform case */
753		− *tx = monolithic_fft;
754		− if (type == AV_TX_FLOAT_MDCT)
755		− *tx = inv ? monolithic_imdct : monolithic_mdct;
756		− }
757		−
758		− if (n != 1)
759		− ff_thread_once(&tabs_53_once, ff_init_53_tabs);
760		− if (m != 1) {
761		− get_ptwo_revtab(s, m, inv);
762		− for (int i = 4; i <= av_log2(m); i++)
763		− ff_init_ff_cos_tabs(i);
764		− }
765		−
766		− if (type == AV_TX_FLOAT_MDCT)
767		− if ((err = gen_mdct_exptab(s, nm, ((float *)scale))))
768		− return err;
769		−
770		− s->n = n;
771		− s->m = m;
772		−
773		− return 0;
774		−}
775		−
776	112	av_cold int av_tx_init(AVTXContext *ctx, av_tx_fn tx, enum AVTXType type,
777	113	int inv, int len, const void *scale, uint64_t flags)
778	114	{
@@ -784,7 +120,12 @@ av_cold int av_tx_init(AVTXContext *ctx, av_tx_fn tx, enum AVTXType type, · 5 unmodified lines
784	120	switch (type) {
785	121	case AV_TX_FLOAT_FFT:
786	122	case AV_TX_FLOAT_MDCT:
787		− if ((err = init_mdct_fft(s, tx, type, inv, len, scale, flags)))
	123	+ if ((err = ff_tx_init_mdct_fft_float(s, tx, type, inv, len, scale, flags)))
	124	+ goto fail;
	125	+ break;
	126	+ case AV_TX_DOUBLE_FFT:
	127	+ case AV_TX_DOUBLE_MDCT:
	128	+ if ((err = ff_tx_init_mdct_fft_double(s, tx, type, inv, len, scale, flags)))
788	129	goto fail;
789	130	break;
790	131	default:

libavutil/tx.h+12−0

@@ -28,6 +28,10 @@ typedef struct AVComplexFloat { · 27 unmodified lines
28	28	float re, im;
29	29	} AVComplexFloat;
30	30
	31	+typedef struct AVComplexDouble {
	32	+ double re, im;
	33	+} AVComplexDouble;
	34	+
31	35	enum AVTXType {
32	36	/**
33	37	* Standard complex to complex FFT with sample data type AVComplexFloat.
@@ -39,6 +43,14 @@ enum AVTXType { · 5 unmodified lines
39	43	* float. Length is the frame size, not the window size (which is 2x frame)
40	44	*/
41	45	AV_TX_FLOAT_MDCT = 1,
	46	+ /**
	47	+ * Same as AV_TX_FLOAT_FFT with a data type of AVComplexDouble.
	48	+ */
	49	+ AV_TX_DOUBLE_FFT = 2,
	50	+ /**
	51	+ * Same as AV_TX_FLOAT_MDCT with data and scale type of double.
	52	+ */
	53	+ AV_TX_DOUBLE_MDCT = 3,
42	54	};
43	55
44	56	/**

libavutil/tx_double.c+21−0

@@ -0,0 +1,21 @@
	1	+/*
	2	+ * This file is part of FFmpeg.
	3	+ *
	4	+ * FFmpeg is free software; you can redistribute it and/or
	5	+ * modify it under the terms of the GNU Lesser General Public
	6	+ * License as published by the Free Software Foundation; either
	7	+ * version 2.1 of the License, or (at your option) any later version.
	8	+ *
	9	+ * FFmpeg is distributed in the hope that it will be useful,
	10	+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
	11	+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	12	+ * Lesser General Public License for more details.
	13	+ *
	14	+ * You should have received a copy of the GNU Lesser General Public
	15	+ * License along with FFmpeg; if not, write to the Free Software
	16	+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
	17	+ */
	18	+
	19	+#define TX_DOUBLE
	20	+#include "tx_priv.h"
	21	+#include "tx_template.c"

libavutil/tx_float.c+21−0

@@ -0,0 +1,21 @@
	1	+/*
	2	+ * This file is part of FFmpeg.
	3	+ *
	4	+ * FFmpeg is free software; you can redistribute it and/or
	5	+ * modify it under the terms of the GNU Lesser General Public
	6	+ * License as published by the Free Software Foundation; either
	7	+ * version 2.1 of the License, or (at your option) any later version.
	8	+ *
	9	+ * FFmpeg is distributed in the hope that it will be useful,
	10	+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
	11	+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	12	+ * Lesser General Public License for more details.
	13	+ *
	14	+ * You should have received a copy of the GNU Lesser General Public
	15	+ * License along with FFmpeg; if not, write to the Free Software
	16	+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
	17	+ */
	18	+
	19	+#define TX_FLOAT
	20	+#include "tx_priv.h"
	21	+#include "tx_template.c"

libavutil/tx_priv.h+105−0

@@ -0,0 +1,105 @@
	1	+/*
	2	+ * This file is part of FFmpeg.
	3	+ *
	4	+ * FFmpeg is free software; you can redistribute it and/or
	5	+ * modify it under the terms of the GNU Lesser General Public
	6	+ * License as published by the Free Software Foundation; either
	7	+ * version 2.1 of the License, or (at your option) any later version.
	8	+ *
	9	+ * FFmpeg is distributed in the hope that it will be useful,
	10	+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
	11	+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	12	+ * Lesser General Public License for more details.
	13	+ *
	14	+ * You should have received a copy of the GNU Lesser General Public
	15	+ * License along with FFmpeg; if not, write to the Free Software
	16	+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
	17	+ */
	18	+
	19	+#ifndef AVUTIL_TX_PRIV_H
	20	+#define AVUTIL_TX_PRIV_H
	21	+
	22	+#include "tx.h"
	23	+#include <stddef.h>
	24	+#include "thread.h"
	25	+#include "mem.h"
	26	+#include "avassert.h"
	27	+#include "attributes.h"
	28	+
	29	+#ifdef TX_FLOAT
	30	+#define TX_NAME(x) x ## _float
	31	+typedef float FFTSample;
	32	+typedef AVComplexFloat FFTComplex;
	33	+#elif defined(TX_DOUBLE)
	34	+#define TX_NAME(x) x ## _double
	35	+typedef double FFTSample;
	36	+typedef AVComplexDouble FFTComplex;
	37	+#else
	38	+typedef void FFTComplex;
	39	+#endif
	40	+
	41	+#if defined(TX_FLOAT) \|\| defined(TX_DOUBLE)
	42	+#define BF(x, y, a, b) do { \
	43	+ x = (a) - (b); \
	44	+ y = (a) + (b); \
	45	+ } while (0)
	46	+
	47	+#define CMUL(dre, dim, are, aim, bre, bim) do { \
	48	+ (dre) = (are) * (bre) - (aim) * (bim); \
	49	+ (dim) = (are) * (bim) + (aim) * (bre); \
	50	+ } while (0)
	51	+#endif
	52	+
	53	+#define CMUL3(c, a, b) \
	54	+ CMUL((c).re, (c).im, (a).re, (a).im, (b).re, (b).im)
	55	+
	56	+#define COSTABLE(size) \
	57	+ DECLARE_ALIGNED(32, FFTSample, TX_NAME(ff_cos_##size))[size/2]
	58	+
	59	+/* Used by asm, reorder with care */
	60	+struct AVTXContext {
	61	+ int n; /* Nptwo part */
	62	+ int m; /* Ptwo part */
	63	+ int inv; /* Is inverted */
	64	+ int type; /* Type */
	65	+
	66	+ FFTComplex exptab; / MDCT exptab */
	67	+ FFTComplex tmp; / Temporary buffer needed for all compound transforms */
	68	+ int pfatab; / Input/Output mapping for compound transforms */
	69	+ int revtab; / Input mapping for power of two transforms */
	70	+};
	71	+
	72	+/* Shared functions */
	73	+int ff_tx_gen_compound_mapping(AVTXContext *s);
	74	+int ff_tx_gen_ptwo_revtab(AVTXContext *s);
	75	+
	76	+/* Also used by SIMD init */
	77	+static inline int split_radix_permutation(int i, int n, int inverse)
	78	+{
	79	+ int m;
	80	+ if (n <= 2)
	81	+ return i & 1;
	82	+ m = n >> 1;
	83	+ if (!(i & m))
	84	+ return split_radix_permutation(i, m, inverse)*2;
	85	+ m >>= 1;
	86	+ if (inverse == !(i & m))
	87	+ return split_radix_permutation(i, m, inverse)*4 + 1;
	88	+ else
	89	+ return split_radix_permutation(i, m, inverse)*4 - 1;
	90	+}
	91	+
	92	+/* Templated functions */
	93	+int ff_tx_init_mdct_fft_float(AVTXContext s, av_tx_fn tx,
	94	+ enum AVTXType type, int inv, int len,
	95	+ const void *scale, uint64_t flags);
	96	+int ff_tx_init_mdct_fft_double(AVTXContext s, av_tx_fn tx,
	97	+ enum AVTXType type, int inv, int len,
	98	+ const void *scale, uint64_t flags);
	99	+
	100	+typedef struct CosTabsInitOnce {
	101	+ void (*func)(void);
	102	+ AVOnce control;
	103	+} CosTabsInitOnce;
	104	+
	105	+#endif /* AVUTIL_TX_PRIV_H */

libavutil/tx_template.c+643−0

@@ -0,0 +1,643 @@
	1	+/*
	2	+ * Copyright (c) 2019 Lynne <dev@lynne.ee>
	3	+ * Power of two FFT:
	4	+ * Copyright (c) 2008 Loren Merritt
	5	+ * Copyright (c) 2002 Fabrice Bellard
	6	+ * Partly based on libdjbfft by D. J. Bernstein
	7	+ *
	8	+ * This file is part of FFmpeg.
	9	+ *
	10	+ * FFmpeg is free software; you can redistribute it and/or
	11	+ * modify it under the terms of the GNU Lesser General Public
	12	+ * License as published by the Free Software Foundation; either
	13	+ * version 2.1 of the License, or (at your option) any later version.
	14	+ *
	15	+ * FFmpeg is distributed in the hope that it will be useful,
	16	+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
	17	+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	18	+ * Lesser General Public License for more details.
	19	+ *
	20	+ * You should have received a copy of the GNU Lesser General Public
	21	+ * License along with FFmpeg; if not, write to the Free Software
	22	+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
	23	+ */
	24	+
	25	+/* All costabs for a type are defined here */
	26	+COSTABLE(16);
	27	+COSTABLE(32);
	28	+COSTABLE(64);
	29	+COSTABLE(128);
	30	+COSTABLE(256);
	31	+COSTABLE(512);
	32	+COSTABLE(1024);
	33	+COSTABLE(2048);
	34	+COSTABLE(4096);
	35	+COSTABLE(8192);
	36	+COSTABLE(16384);
	37	+COSTABLE(32768);
	38	+COSTABLE(65536);
	39	+COSTABLE(131072);
	40	+DECLARE_ALIGNED(32, FFTComplex, TX_NAME(ff_cos_53))[4];
	41	+
	42	+static FFTSample * const cos_tabs[18] = {
	43	+ NULL,
	44	+ NULL,
	45	+ NULL,
	46	+ NULL,
	47	+ TX_NAME(ff_cos_16),
	48	+ TX_NAME(ff_cos_32),
	49	+ TX_NAME(ff_cos_64),
	50	+ TX_NAME(ff_cos_128),
	51	+ TX_NAME(ff_cos_256),
	52	+ TX_NAME(ff_cos_512),
	53	+ TX_NAME(ff_cos_1024),
	54	+ TX_NAME(ff_cos_2048),
	55	+ TX_NAME(ff_cos_4096),
	56	+ TX_NAME(ff_cos_8192),
	57	+ TX_NAME(ff_cos_16384),
	58	+ TX_NAME(ff_cos_32768),
	59	+ TX_NAME(ff_cos_65536),
	60	+ TX_NAME(ff_cos_131072),
	61	+};
	62	+
	63	+static av_always_inline void init_cos_tabs_idx(int index)
	64	+{
	65	+ int m = 1 << index;
	66	+ double freq = 2*M_PI/m;
	67	+ FFTSample *tab = cos_tabs[index];
	68	+ for(int i = 0; i <= m/4; i++)
	69	+ tab[i] = cos(i*freq);
	70	+ for(int i = 1; i < m/4; i++)
	71	+ tab[m/2 - i] = tab[i];
	72	+}
	73	+
	74	+#define INIT_FF_COS_TABS_FUNC(index, size) \
	75	+static av_cold void init_cos_tabs_ ## size (void) \
	76	+{ \
	77	+ init_cos_tabs_idx(index); \
	78	+}
	79	+
	80	+INIT_FF_COS_TABS_FUNC(4, 16)
	81	+INIT_FF_COS_TABS_FUNC(5, 32)
	82	+INIT_FF_COS_TABS_FUNC(6, 64)
	83	+INIT_FF_COS_TABS_FUNC(7, 128)
	84	+INIT_FF_COS_TABS_FUNC(8, 256)
	85	+INIT_FF_COS_TABS_FUNC(9, 512)
	86	+INIT_FF_COS_TABS_FUNC(10, 1024)
	87	+INIT_FF_COS_TABS_FUNC(11, 2048)
	88	+INIT_FF_COS_TABS_FUNC(12, 4096)
	89	+INIT_FF_COS_TABS_FUNC(13, 8192)
	90	+INIT_FF_COS_TABS_FUNC(14, 16384)
	91	+INIT_FF_COS_TABS_FUNC(15, 32768)
	92	+INIT_FF_COS_TABS_FUNC(16, 65536)
	93	+INIT_FF_COS_TABS_FUNC(17, 131072)
	94	+
	95	+static av_cold void ff_init_53_tabs(void)
	96	+{
	97	+ TX_NAME(ff_cos_53)[0] = (FFTComplex){ cos(2 * M_PI / 12), cos(2 * M_PI / 12) };
	98	+ TX_NAME(ff_cos_53)[1] = (FFTComplex){ 0.5, 0.5 };
	99	+ TX_NAME(ff_cos_53)[2] = (FFTComplex){ cos(2 * M_PI / 5), sin(2 * M_PI / 5) };
	100	+ TX_NAME(ff_cos_53)[3] = (FFTComplex){ cos(2 * M_PI / 10), sin(2 * M_PI / 10) };
	101	+}
	102	+
	103	+static CosTabsInitOnce cos_tabs_init_once[] = {
	104	+ { ff_init_53_tabs, AV_ONCE_INIT },
	105	+ { NULL },
	106	+ { NULL },
	107	+ { NULL },
	108	+ { init_cos_tabs_16, AV_ONCE_INIT },
	109	+ { init_cos_tabs_32, AV_ONCE_INIT },
	110	+ { init_cos_tabs_64, AV_ONCE_INIT },
	111	+ { init_cos_tabs_128, AV_ONCE_INIT },
	112	+ { init_cos_tabs_256, AV_ONCE_INIT },
	113	+ { init_cos_tabs_512, AV_ONCE_INIT },
	114	+ { init_cos_tabs_1024, AV_ONCE_INIT },
	115	+ { init_cos_tabs_2048, AV_ONCE_INIT },
	116	+ { init_cos_tabs_4096, AV_ONCE_INIT },
	117	+ { init_cos_tabs_8192, AV_ONCE_INIT },
	118	+ { init_cos_tabs_16384, AV_ONCE_INIT },
	119	+ { init_cos_tabs_32768, AV_ONCE_INIT },
	120	+ { init_cos_tabs_65536, AV_ONCE_INIT },
	121	+ { init_cos_tabs_131072, AV_ONCE_INIT },
	122	+};
	123	+
	124	+static av_cold void init_cos_tabs(int index)
	125	+{
	126	+ ff_thread_once(&cos_tabs_init_once[index].control,
	127	+ cos_tabs_init_once[index].func);
	128	+}
	129	+
	130	+static av_always_inline void fft3(FFTComplex out, FFTComplex in,
	131	+ ptrdiff_t stride)
	132	+{
	133	+ FFTComplex tmp[2];
	134	+
	135	+ tmp[0].re = in[1].im - in[2].im;
	136	+ tmp[0].im = in[1].re - in[2].re;
	137	+ tmp[1].re = in[1].re + in[2].re;
	138	+ tmp[1].im = in[1].im + in[2].im;
	139	+
	140	+ out[0*stride].re = in[0].re + tmp[1].re;
	141	+ out[0*stride].im = in[0].im + tmp[1].im;
	142	+
	143	+ tmp[0].re *= TX_NAME(ff_cos_53)[0].re;
	144	+ tmp[0].im *= TX_NAME(ff_cos_53)[0].im;
	145	+ tmp[1].re *= TX_NAME(ff_cos_53)[1].re;
	146	+ tmp[1].im *= TX_NAME(ff_cos_53)[1].re;
	147	+
	148	+ out[1*stride].re = in[0].re - tmp[1].re + tmp[0].re;
	149	+ out[1*stride].im = in[0].im - tmp[1].im - tmp[0].im;
	150	+ out[2*stride].re = in[0].re - tmp[1].re - tmp[0].re;
	151	+ out[2*stride].im = in[0].im - tmp[1].im + tmp[0].im;
	152	+}
	153	+
	154	+#define DECL_FFT5(NAME, D0, D1, D2, D3, D4) \
	155	+static av_always_inline void NAME(FFTComplex out, FFTComplex in, \
	156	+ ptrdiff_t stride) \
	157	+{ \
	158	+ FFTComplex z0[4], t[6]; \
	159	+ \
	160	+ t[0].re = in[1].re + in[4].re; \
	161	+ t[0].im = in[1].im + in[4].im; \
	162	+ t[1].im = in[1].re - in[4].re; \
	163	+ t[1].re = in[1].im - in[4].im; \
	164	+ t[2].re = in[2].re + in[3].re; \
	165	+ t[2].im = in[2].im + in[3].im; \
	166	+ t[3].im = in[2].re - in[3].re; \
	167	+ t[3].re = in[2].im - in[3].im; \
	168	+ \
	169	+ out[D0*stride].re = in[0].re + in[1].re + in[2].re + \
	170	+ in[3].re + in[4].re; \
	171	+ out[D0*stride].im = in[0].im + in[1].im + in[2].im + \
	172	+ in[3].im + in[4].im; \
	173	+ \
	174	+ t[4].re = TX_NAME(ff_cos_53)[2].re * t[2].re; \
	175	+ t[4].im = TX_NAME(ff_cos_53)[2].re * t[2].im; \
	176	+ t[4].re -= TX_NAME(ff_cos_53)[3].re * t[0].re; \
	177	+ t[4].im -= TX_NAME(ff_cos_53)[3].re * t[0].im; \
	178	+ t[0].re = TX_NAME(ff_cos_53)[2].re * t[0].re; \
	179	+ t[0].im = TX_NAME(ff_cos_53)[2].re * t[0].im; \
	180	+ t[0].re -= TX_NAME(ff_cos_53)[3].re * t[2].re; \
	181	+ t[0].im -= TX_NAME(ff_cos_53)[3].re * t[2].im; \
	182	+ t[5].re = TX_NAME(ff_cos_53)[2].im * t[3].re; \
	183	+ t[5].im = TX_NAME(ff_cos_53)[2].im * t[3].im; \
	184	+ t[5].re -= TX_NAME(ff_cos_53)[3].im * t[1].re; \
	185	+ t[5].im -= TX_NAME(ff_cos_53)[3].im * t[1].im; \
	186	+ t[1].re = TX_NAME(ff_cos_53)[2].im * t[1].re; \
	187	+ t[1].im = TX_NAME(ff_cos_53)[2].im * t[1].im; \
	188	+ t[1].re += TX_NAME(ff_cos_53)[3].im * t[3].re; \
	189	+ t[1].im += TX_NAME(ff_cos_53)[3].im * t[3].im; \
	190	+ \
	191	+ z0[0].re = t[0].re - t[1].re; \
	192	+ z0[0].im = t[0].im - t[1].im; \
	193	+ z0[1].re = t[4].re + t[5].re; \
	194	+ z0[1].im = t[4].im + t[5].im; \
	195	+ \
	196	+ z0[2].re = t[4].re - t[5].re; \
	197	+ z0[2].im = t[4].im - t[5].im; \
	198	+ z0[3].re = t[0].re + t[1].re; \
	199	+ z0[3].im = t[0].im + t[1].im; \
	200	+ \
	201	+ out[D1*stride].re = in[0].re + z0[3].re; \
	202	+ out[D1*stride].im = in[0].im + z0[0].im; \
	203	+ out[D2*stride].re = in[0].re + z0[2].re; \
	204	+ out[D2*stride].im = in[0].im + z0[1].im; \
	205	+ out[D3*stride].re = in[0].re + z0[1].re; \
	206	+ out[D3*stride].im = in[0].im + z0[2].im; \
	207	+ out[D4*stride].re = in[0].re + z0[0].re; \
	208	+ out[D4*stride].im = in[0].im + z0[3].im; \
	209	+}
	210	+
	211	+DECL_FFT5(fft5, 0, 1, 2, 3, 4)
	212	+DECL_FFT5(fft5_m1, 0, 6, 12, 3, 9)
	213	+DECL_FFT5(fft5_m2, 10, 1, 7, 13, 4)
	214	+DECL_FFT5(fft5_m3, 5, 11, 2, 8, 14)
	215	+
	216	+static av_always_inline void fft15(FFTComplex out, FFTComplex in,
	217	+ ptrdiff_t stride)
	218	+{
	219	+ FFTComplex tmp[15];
	220	+
	221	+ for (int i = 0; i < 5; i++)
	222	+ fft3(tmp + i, in + i*3, 5);
	223	+
	224	+ fft5_m1(out, tmp + 0, stride);
	225	+ fft5_m2(out, tmp + 5, stride);
	226	+ fft5_m3(out, tmp + 10, stride);
	227	+}
	228	+
	229	+#define BUTTERFLIES(a0,a1,a2,a3) {\
	230	+ BF(t3, t5, t5, t1);\
	231	+ BF(a2.re, a0.re, a0.re, t5);\
	232	+ BF(a3.im, a1.im, a1.im, t3);\
	233	+ BF(t4, t6, t2, t6);\
	234	+ BF(a3.re, a1.re, a1.re, t4);\
	235	+ BF(a2.im, a0.im, a0.im, t6);\
	236	+}
	237	+
	238	+// force loading all the inputs before storing any.
	239	+// this is slightly slower for small data, but avoids store->load aliasing
	240	+// for addresses separated by large powers of 2.
	241	+#define BUTTERFLIES_BIG(a0,a1,a2,a3) {\
	242	+ FFTSample r0=a0.re, i0=a0.im, r1=a1.re, i1=a1.im;\
	243	+ BF(t3, t5, t5, t1);\
	244	+ BF(a2.re, a0.re, r0, t5);\
	245	+ BF(a3.im, a1.im, i1, t3);\
	246	+ BF(t4, t6, t2, t6);\
	247	+ BF(a3.re, a1.re, r1, t4);\
	248	+ BF(a2.im, a0.im, i0, t6);\
	249	+}
	250	+
	251	+#define TRANSFORM(a0,a1,a2,a3,wre,wim) {\
	252	+ CMUL(t1, t2, a2.re, a2.im, wre, -wim);\
	253	+ CMUL(t5, t6, a3.re, a3.im, wre, wim);\
	254	+ BUTTERFLIES(a0,a1,a2,a3)\
	255	+}
	256	+
	257	+#define TRANSFORM_ZERO(a0,a1,a2,a3) {\
	258	+ t1 = a2.re;\
	259	+ t2 = a2.im;\
	260	+ t5 = a3.re;\
	261	+ t6 = a3.im;\
	262	+ BUTTERFLIES(a0,a1,a2,a3)\
	263	+}
	264	+
	265	+/* z[0...8n-1], w[1...2n-1] */
	266	+#define PASS(name)\
	267	+static void name(FFTComplex z, const FFTSample wre, unsigned int n)\
	268	+{\
	269	+ FFTSample t1, t2, t3, t4, t5, t6;\
	270	+ int o1 = 2*n;\
	271	+ int o2 = 4*n;\
	272	+ int o3 = 6*n;\
	273	+ const FFTSample *wim = wre+o1;\
	274	+ n--;\
	275	+\
	276	+ TRANSFORM_ZERO(z[0],z[o1],z[o2],z[o3]);\
	277	+ TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
	278	+ do {\
	279	+ z += 2;\
	280	+ wre += 2;\
	281	+ wim -= 2;\
	282	+ TRANSFORM(z[0],z[o1],z[o2],z[o3],wre[0],wim[0]);\
	283	+ TRANSFORM(z[1],z[o1+1],z[o2+1],z[o3+1],wre[1],wim[-1]);\
	284	+ } while(--n);\
	285	+}
	286	+
	287	+PASS(pass)
	288	+#undef BUTTERFLIES
	289	+#define BUTTERFLIES BUTTERFLIES_BIG
	290	+PASS(pass_big)
	291	+
	292	+#define DECL_FFT(n,n2,n4)\
	293	+static void fft##n(FFTComplex *z)\
	294	+{\
	295	+ fft##n2(z);\
	296	+ fft##n4(z+n4*2);\
	297	+ fft##n4(z+n4*3);\
	298	+ pass(z,TX_NAME(ff_cos_##n),n4/2);\
	299	+}
	300	+
	301	+static void fft4(FFTComplex *z)
	302	+{
	303	+ FFTSample t1, t2, t3, t4, t5, t6, t7, t8;
	304	+
	305	+ BF(t3, t1, z[0].re, z[1].re);
	306	+ BF(t8, t6, z[3].re, z[2].re);
	307	+ BF(z[2].re, z[0].re, t1, t6);
	308	+ BF(t4, t2, z[0].im, z[1].im);
	309	+ BF(t7, t5, z[2].im, z[3].im);
	310	+ BF(z[3].im, z[1].im, t4, t8);
	311	+ BF(z[3].re, z[1].re, t3, t7);
	312	+ BF(z[2].im, z[0].im, t2, t5);
	313	+}
	314	+
	315	+static void fft8(FFTComplex *z)
	316	+{
	317	+ FFTSample t1, t2, t3, t4, t5, t6;
	318	+
	319	+ fft4(z);
	320	+
	321	+ BF(t1, z[5].re, z[4].re, -z[5].re);
	322	+ BF(t2, z[5].im, z[4].im, -z[5].im);
	323	+ BF(t5, z[7].re, z[6].re, -z[7].re);
	324	+ BF(t6, z[7].im, z[6].im, -z[7].im);
	325	+
	326	+ BUTTERFLIES(z[0],z[2],z[4],z[6]);
	327	+ TRANSFORM(z[1],z[3],z[5],z[7],M_SQRT1_2,M_SQRT1_2);
	328	+}
	329	+
	330	+static void fft16(FFTComplex *z)
	331	+{
	332	+ FFTSample t1, t2, t3, t4, t5, t6;
	333	+ FFTSample cos_16_1 = TX_NAME(ff_cos_16)[1];
	334	+ FFTSample cos_16_3 = TX_NAME(ff_cos_16)[3];
	335	+
	336	+ fft8(z);
	337	+ fft4(z+8);
	338	+ fft4(z+12);
	339	+
	340	+ TRANSFORM_ZERO(z[0],z[4],z[8],z[12]);
	341	+ TRANSFORM(z[2],z[6],z[10],z[14],M_SQRT1_2,M_SQRT1_2);
	342	+ TRANSFORM(z[1],z[5],z[9],z[13],cos_16_1,cos_16_3);
	343	+ TRANSFORM(z[3],z[7],z[11],z[15],cos_16_3,cos_16_1);
	344	+}
	345	+
	346	+DECL_FFT(32,16,8)
	347	+DECL_FFT(64,32,16)
	348	+DECL_FFT(128,64,32)
	349	+DECL_FFT(256,128,64)
	350	+DECL_FFT(512,256,128)
	351	+#define pass pass_big
	352	+DECL_FFT(1024,512,256)
	353	+DECL_FFT(2048,1024,512)
	354	+DECL_FFT(4096,2048,1024)
	355	+DECL_FFT(8192,4096,2048)
	356	+DECL_FFT(16384,8192,4096)
	357	+DECL_FFT(32768,16384,8192)
	358	+DECL_FFT(65536,32768,16384)
	359	+DECL_FFT(131072,65536,32768)
	360	+
	361	+static void (* const fft_dispatch[])(FFTComplex*) = {
	362	+ fft4, fft8, fft16, fft32, fft64, fft128, fft256, fft512, fft1024,
	363	+ fft2048, fft4096, fft8192, fft16384, fft32768, fft65536, fft131072
	364	+};
	365	+
	366	+#define DECL_COMP_FFT(N) \
	367	+static void compound_fft_##N##xM(AVTXContext s, void _out, \
	368	+ void *_in, ptrdiff_t stride) \
	369	+{ \
	370	+ const int m = s->m, in_map = s->pfatab, out_map = in_map + N*m; \
	371	+ FFTComplex *in = _in; \
	372	+ FFTComplex *out = _out; \
	373	+ FFTComplex fft##N##in[N]; \
	374	+ void (fftp)(FFTComplex z) = fft_dispatch[av_log2(m) - 2]; \
	375	+ \
	376	+ for (int i = 0; i < m; i++) { \
	377	+ for (int j = 0; j < N; j++) \
	378	+ fft##N##in[j] = in[in_map[i*N + j]]; \
	379	+ fft##N(s->tmp + s->revtab[i], fft##N##in, m); \
	380	+ } \
	381	+ \
	382	+ for (int i = 0; i < N; i++) \
	383	+ fftp(s->tmp + m*i); \
	384	+ \
	385	+ for (int i = 0; i < N*m; i++) \
	386	+ out[i] = s->tmp[out_map[i]]; \
	387	+}
	388	+
	389	+DECL_COMP_FFT(3)
	390	+DECL_COMP_FFT(5)
	391	+DECL_COMP_FFT(15)
	392	+
	393	+static void monolithic_fft(AVTXContext s, void _out, void *_in,
	394	+ ptrdiff_t stride)
	395	+{
	396	+ FFTComplex *in = _in;
	397	+ FFTComplex *out = _out;
	398	+ int m = s->m, mb = av_log2(m) - 2;
	399	+ for (int i = 0; i < m; i++)
	400	+ out[s->revtab[i]] = in[i];
	401	+ fft_dispatch[mb](out);
	402	+}
	403	+
	404	+#define DECL_COMP_IMDCT(N) \
	405	+static void compound_imdct_##N##xM(AVTXContext s, void _dst, void *_src, \
	406	+ ptrdiff_t stride) \
	407	+{ \
	408	+ FFTComplex fft##N##in[N]; \
	409	+ FFTComplex z = _dst, exp = s->exptab; \
	410	+ const int m = s->m, len8 = N*m >> 1; \
	411	+ const int in_map = s->pfatab, out_map = in_map + N*m; \
	412	+ const FFTSample src = _src, in1, *in2; \
	413	+ void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2]; \
	414	+ \
	415	+ stride /= sizeof(src); / To convert it from bytes */ \
	416	+ in1 = src; \
	417	+ in2 = src + ((Nm2) - 1) * stride; \
	418	+ \
	419	+ for (int i = 0; i < m; i++) { \
	420	+ for (int j = 0; j < N; j++) { \
	421	+ const int k = in_map[i*N + j]; \
	422	+ FFTComplex tmp = { in2[-kstride], in1[kstride] }; \
	423	+ CMUL3(fft##N##in[j], tmp, exp[k >> 1]); \
	424	+ } \
	425	+ fft##N(s->tmp + s->revtab[i], fft##N##in, m); \
	426	+ } \
	427	+ \
	428	+ for (int i = 0; i < N; i++) \
	429	+ fftp(s->tmp + m*i); \
	430	+ \
	431	+ for (int i = 0; i < len8; i++) { \
	432	+ const int i0 = len8 + i, i1 = len8 - i - 1; \
	433	+ const int s0 = out_map[i0], s1 = out_map[i1]; \
	434	+ FFTComplex src1 = { s->tmp[s1].im, s->tmp[s1].re }; \
	435	+ FFTComplex src0 = { s->tmp[s0].im, s->tmp[s0].re }; \
	436	+ \
	437	+ CMUL(z[i1].re, z[i0].im, src1.re, src1.im, exp[i1].im, exp[i1].re); \
	438	+ CMUL(z[i0].re, z[i1].im, src0.re, src0.im, exp[i0].im, exp[i0].re); \
	439	+ } \
	440	+}
	441	+
	442	+DECL_COMP_IMDCT(3)
	443	+DECL_COMP_IMDCT(5)
	444	+DECL_COMP_IMDCT(15)
	445	+
	446	+#define DECL_COMP_MDCT(N) \
	447	+static void compound_mdct_##N##xM(AVTXContext s, void _dst, void *_src, \
	448	+ ptrdiff_t stride) \
	449	+{ \
	450	+ FFTSample src = _src, dst = _dst; \
	451	+ FFTComplex *exp = s->exptab, tmp, fft##N##in[N]; \
	452	+ const int m = s->m, len4 = Nm, len3 = len4 3, len8 = len4 >> 1; \
	453	+ const int in_map = s->pfatab, out_map = in_map + N*m; \
	454	+ void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2]; \
	455	+ \
	456	+ stride /= sizeof(*dst); \
	457	+ \
	458	+ for (int i = 0; i < m; i++) { /* Folding and pre-reindexing */ \
	459	+ for (int j = 0; j < N; j++) { \
	460	+ const int k = in_map[i*N + j]; \
	461	+ if (k < len4) { \
	462	+ tmp.re = -src[ len4 + k] + src[1*len4 - 1 - k]; \
	463	+ tmp.im = -src[ len3 + k] - src[1*len3 - 1 - k]; \
	464	+ } else { \
	465	+ tmp.re = -src[ len4 + k] - src[5*len4 - 1 - k]; \
	466	+ tmp.im = src[-len4 + k] - src[1*len3 - 1 - k]; \
	467	+ } \
	468	+ CMUL(fft##N##in[j].im, fft##N##in[j].re, tmp.re, tmp.im, \
	469	+ exp[k >> 1].re, exp[k >> 1].im); \
	470	+ } \
	471	+ fft##N(s->tmp + s->revtab[i], fft##N##in, m); \
	472	+ } \
	473	+ \
	474	+ for (int i = 0; i < N; i++) \
	475	+ fftp(s->tmp + m*i); \
	476	+ \
	477	+ for (int i = 0; i < len8; i++) { \
	478	+ const int i0 = len8 + i, i1 = len8 - i - 1; \
	479	+ const int s0 = out_map[i0], s1 = out_map[i1]; \
	480	+ FFTComplex src1 = { s->tmp[s1].re, s->tmp[s1].im }; \
	481	+ FFTComplex src0 = { s->tmp[s0].re, s->tmp[s0].im }; \
	482	+ \
	483	+ CMUL(dst[2i1stride + stride], dst[2i0stride], src0.re, src0.im, \
	484	+ exp[i0].im, exp[i0].re); \
	485	+ CMUL(dst[2i0stride + stride], dst[2i1stride], src1.re, src1.im, \
	486	+ exp[i1].im, exp[i1].re); \
	487	+ } \
	488	+}
	489	+
	490	+DECL_COMP_MDCT(3)
	491	+DECL_COMP_MDCT(5)
	492	+DECL_COMP_MDCT(15)
	493	+
	494	+static void monolithic_imdct(AVTXContext s, void _dst, void *_src,
	495	+ ptrdiff_t stride)
	496	+{
	497	+ FFTComplex z = _dst, exp = s->exptab;
	498	+ const int m = s->m, len8 = m >> 1;
	499	+ const FFTSample src = _src, in1, *in2;
	500	+ void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2];
	501	+
	502	+ stride /= sizeof(*src);
	503	+ in1 = src;
	504	+ in2 = src + ((m2) - 1) stride;
	505	+
	506	+ for (int i = 0; i < m; i++) {
	507	+ FFTComplex tmp = { in2[-2istride], in1[2istride] };
	508	+ CMUL3(z[s->revtab[i]], tmp, exp[i]);
	509	+ }
	510	+
	511	+ fftp(z);
	512	+
	513	+ for (int i = 0; i < len8; i++) {
	514	+ const int i0 = len8 + i, i1 = len8 - i - 1;
	515	+ FFTComplex src1 = { z[i1].im, z[i1].re };
	516	+ FFTComplex src0 = { z[i0].im, z[i0].re };
	517	+
	518	+ CMUL(z[i1].re, z[i0].im, src1.re, src1.im, exp[i1].im, exp[i1].re);
	519	+ CMUL(z[i0].re, z[i1].im, src0.re, src0.im, exp[i0].im, exp[i0].re);
	520	+ }
	521	+}
	522	+
	523	+static void monolithic_mdct(AVTXContext s, void _dst, void *_src,
	524	+ ptrdiff_t stride)
	525	+{
	526	+ FFTSample src = _src, dst = _dst;
	527	+ FFTComplex exp = s->exptab, tmp, z = _dst;
	528	+ const int m = s->m, len4 = m, len3 = len4 * 3, len8 = len4 >> 1;
	529	+ void (fftp)(FFTComplex ) = fft_dispatch[av_log2(m) - 2];
	530	+
	531	+ stride /= sizeof(*dst);
	532	+
	533	+ for (int i = 0; i < m; i++) { /* Folding and pre-reindexing */
	534	+ const int k = 2*i;
	535	+ if (k < len4) {
	536	+ tmp.re = -src[ len4 + k] + src[1*len4 - 1 - k];
	537	+ tmp.im = -src[ len3 + k] - src[1*len3 - 1 - k];
	538	+ } else {
	539	+ tmp.re = -src[ len4 + k] - src[5*len4 - 1 - k];
	540	+ tmp.im = src[-len4 + k] - src[1*len3 - 1 - k];
	541	+ }
	542	+ CMUL(z[s->revtab[i]].im, z[s->revtab[i]].re, tmp.re, tmp.im,
	543	+ exp[i].re, exp[i].im);
	544	+ }
	545	+
	546	+ fftp(z);
	547	+
	548	+ for (int i = 0; i < len8; i++) {
	549	+ const int i0 = len8 + i, i1 = len8 - i - 1;
	550	+ FFTComplex src1 = { z[i1].re, z[i1].im };
	551	+ FFTComplex src0 = { z[i0].re, z[i0].im };
	552	+
	553	+ CMUL(dst[2i1stride + stride], dst[2i0stride], src0.re, src0.im,
	554	+ exp[i0].im, exp[i0].re);
	555	+ CMUL(dst[2i0stride + stride], dst[2i1stride], src1.re, src1.im,
	556	+ exp[i1].im, exp[i1].re);
	557	+ }
	558	+}
	559	+
	560	+static int gen_mdct_exptab(AVTXContext *s, int len4, double scale)
	561	+{
	562	+ const double theta = (scale < 0 ? len4 : 0) + 1.0/8.0;
	563	+
	564	+ if (!(s->exptab = av_malloc_array(len4, sizeof(*s->exptab))))
	565	+ return AVERROR(ENOMEM);
	566	+
	567	+ scale = sqrt(fabs(scale));
	568	+ for (int i = 0; i < len4; i++) {
	569	+ const double alpha = M_PI_2 * (i + theta) / len4;
	570	+ s->exptab[i].re = cos(alpha) * scale;
	571	+ s->exptab[i].im = sin(alpha) * scale;
	572	+ }
	573	+
	574	+ return 0;
	575	+}
	576	+
	577	+int TX_NAME(ff_tx_init_mdct_fft)(AVTXContext s, av_tx_fn tx,
	578	+ enum AVTXType type, int inv, int len,
	579	+ const void *scale, uint64_t flags)
	580	+{
	581	+ const int is_mdct = type == AV_TX_FLOAT_MDCT \|\| type == AV_TX_DOUBLE_MDCT;
	582	+ int err, n = 1, m = 1, max_ptwo = 1 << (FF_ARRAY_ELEMS(fft_dispatch) + 1);
	583	+
	584	+ if (is_mdct)
	585	+ len >>= 1;
	586	+
	587	+#define CHECK_FACTOR(DST, FACTOR, SRC) \
	588	+ if (DST == 1 && !(SRC % FACTOR)) { \
	589	+ DST = FACTOR; \
	590	+ SRC /= FACTOR; \
	591	+ }
	592	+ CHECK_FACTOR(n, 15, len)
	593	+ CHECK_FACTOR(n, 5, len)
	594	+ CHECK_FACTOR(n, 3, len)
	595	+#undef CHECK_NPTWO_FACTOR
	596	+
	597	+ /* len must be a power of two now */
	598	+ if (!(len & (len - 1)) && len >= 4 && len <= max_ptwo) {
	599	+ m = len;
	600	+ len = 1;
	601	+ }
	602	+
	603	+ s->n = n;
	604	+ s->m = m;
	605	+ s->inv = inv;
	606	+ s->type = type;
	607	+
	608	+ /* Filter out direct 3, 5 and 15 transforms, too niche */
	609	+ if (len > 1 \|\| m == 1) {
	610	+ av_log(NULL, AV_LOG_ERROR, "Unsupported transform size: n = %i, "
	611	+ "m = %i, residual = %i!\n", n, m, len);
	612	+ return AVERROR(EINVAL);
	613	+ } else if (n > 1 && m > 1) { /* 2D transform case */
	614	+ if ((err = ff_tx_gen_compound_mapping(s)))
	615	+ return err;
	616	+ if (!(s->tmp = av_malloc(nmsizeof(*s->tmp))))
	617	+ return AVERROR(ENOMEM);
	618	+ *tx = n == 3 ? compound_fft_3xM :
	619	+ n == 5 ? compound_fft_5xM :
	620	+ compound_fft_15xM;
	621	+ if (is_mdct)
	622	+ *tx = n == 3 ? inv ? compound_imdct_3xM : compound_mdct_3xM :
	623	+ n == 5 ? inv ? compound_imdct_5xM : compound_mdct_5xM :
	624	+ inv ? compound_imdct_15xM : compound_mdct_15xM;
	625	+ } else { /* Direct transform case */
	626	+ *tx = monolithic_fft;
	627	+ if (is_mdct)
	628	+ *tx = inv ? monolithic_imdct : monolithic_mdct;
	629	+ }
	630	+
	631	+ if (n != 1)
	632	+ init_cos_tabs(0);
	633	+ if (m != 1) {
	634	+ ff_tx_gen_ptwo_revtab(s);
	635	+ for (int i = 4; i <= av_log2(m); i++)
	636	+ init_cos_tabs(i);
	637	+ }
	638	+
	639	+ if (is_mdct)
	640	+ return gen_mdct_exptab(s, nm, ((FFTSample *)scale));
	641	+
	642	+ return 0;
	643	+}

libavutil/version.h+1−1

@@ -79,7 +79,7 @@ · 78 unmodified lines
79	79	*/
80	80
81	81	#define LIBAVUTIL_VERSION_MAJOR 56
82		−#define LIBAVUTIL_VERSION_MINOR 32
	82	+#define LIBAVUTIL_VERSION_MINOR 33
83	83	#define LIBAVUTIL_VERSION_MICRO 100
84	84
85	85	#define LIBAVUTIL_VERSION_INT AV_VERSION_INT(LIBAVUTIL_VERSION_MAJOR, \

lavu/tx: add support for double precision FFT and MDCT

Keyboard shortcuts

General

Issues and pull requests

Issues, pull requests and Inbox

Repository code