x86/lpc: implement a new Welch windowing function

Old one was written with the assumption only even inputs would be given. This very messy replacement supports even and odd inputs, and supports AVX2 for extra speed. The buffers given are usually quite big (4k samples), so the speedup is worth it. The new SSE version is still faster than the old inline asm version by 33%. Also checkasm is provided to make sure this monstrosity works. This fixes some FATE tests.

Lynne committed 4 years ago

Commit 3ade6a81 parent cc367a97 files · +343 −58Browse files

Unified · Split

7 files

libavcodec/x86/Makefile +2 −1
libavcodec/x86/lpc.asm +241 −0
libavcodec/x86/lpc_init.c +15 −57
tests/checkasm/Makefile +1 −0
tests/checkasm/checkasm.c +3 −0
tests/checkasm/checkasm.h +1 −0
tests/checkasm/lpc.c +80 −0

libavcodec/x86/Makefile+2−1

@@ -23,7 +23,7 @@ OBJS-$(CONFIG_LLVIDENCDSP) += x86/lossless_videoencdsp_init.o · 22 unmodified lines
23	23	OBJS-$(CONFIG_HUFFYUVDSP) += x86/huffyuvdsp_init.o
24	24	OBJS-$(CONFIG_HUFFYUVENCDSP) += x86/huffyuvencdsp_init.o
25	25	OBJS-$(CONFIG_IDCTDSP) += x86/idctdsp_init.o
26		−OBJS-$(CONFIG_LPC) += x86/lpc.o
	26	+OBJS-$(CONFIG_LPC) += x86/lpc_init.o
27	27	OBJS-$(CONFIG_MDCT15) += x86/mdct15_init.o
28	28	OBJS-$(CONFIG_ME_CMP) += x86/me_cmp_init.o
29	29	OBJS-$(CONFIG_MPEGAUDIODSP) += x86/mpegaudiodsp.o
@@ -127,6 +127,7 @@ X86ASM-OBJS-$(CONFIG_IDCTDSP) += x86/idctdsp.o · 97 unmodified lines
127	127	X86ASM-OBJS-$(CONFIG_LLAUDDSP) += x86/lossless_audiodsp.o
128	128	X86ASM-OBJS-$(CONFIG_LLVIDDSP) += x86/lossless_videodsp.o
129	129	X86ASM-OBJS-$(CONFIG_LLVIDENCDSP) += x86/lossless_videoencdsp.o
	130	+X86ASM-OBJS-$(CONFIG_LPC) += x86/lpc.o
130	131	X86ASM-OBJS-$(CONFIG_MDCT15) += x86/mdct15.o
131	132	X86ASM-OBJS-$(CONFIG_ME_CMP) += x86/me_cmp.o
132	133	X86ASM-OBJS-$(CONFIG_MPEGAUDIODSP) += x86/imdct36.o

libavcodec/x86/lpc.asm+241−0

@@ -0,0 +1,241 @@
	1	+;******************************************************************************
	2	+;* Copyright (c) Lynne
	3	+;*
	4	+;* This file is part of FFmpeg.
	5	+;*
	6	+;* FFmpeg is free software; you can redistribute it and/or
	7	+;* modify it under the terms of the GNU Lesser General Public
	8	+;* License as published by the Free Software Foundation; either
	9	+;* version 2.1 of the License, or (at your option) any later version.
	10	+;*
	11	+;* FFmpeg is distributed in the hope that it will be useful,
	12	+;* but WITHOUT ANY WARRANTY; without even the implied warranty of
	13	+;* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
	14	+;* Lesser General Public License for more details.
	15	+;*
	16	+;* You should have received a copy of the GNU Lesser General Public
	17	+;* License along with FFmpeg; if not, write to the Free Software
	18	+;* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
	19	+;******************************************************************************
	20	+
	21	+%include "libavutil/x86/x86util.asm"
	22	+
	23	+SECTION_RODATA 32
	24	+
	25	+one_tab: times 4 dq 1.0
	26	+seq_tab_avx2: dq 3.0, 2.0, 1.0, 0.0
	27	+sub_tab: dq -1.0, -2.0, -3.0, -4.0
	28	+add_tab_avx2: times 4 dq 4.0
	29	+dec_tab_avx2: times 4 dq -4.0
	30	+add_tab_sse2: times 2 dq 2.0
	31	+dec_tab_sse2: times 2 dq -2.0
	32	+dec_tab_scalar: times 2 dq -1.0
	33	+seq_tab_sse2: dq 1.0, 0.0
	34	+
	35	+SECTION .text
	36	+
	37	+%macro APPLY_WELCH_FN 0
	38	+cglobal lpc_apply_welch_window, 3, 5, 8, data, len, out, off1, off2
	39	+ cmp lenq, 0
	40	+ je .end
	41	+ cmp lenq, 1
	42	+ je .one
	43	+
	44	+ movapd m6, [one_tab]
	45	+
	46	+ movd xm1, lend
	47	+ cvtdq2pd xm1, xm1 ; len
	48	+%if cpuflag(avx2)
	49	+ vbroadcastsd m1, xm1
	50	+%else
	51	+ shufpd m1, m1, 00b
	52	+%endif
	53	+
	54	+ addpd m0, m6, m6 ; 2.0
	55	+ subpd m1, m6 ; len - 1
	56	+ divpd m0, m1 ; 2.0 / (len - 1)
	57	+
	58	+ mov off1q, lenq
	59	+ and off1q, 1
	60	+ je .even
	61	+
	62	+ movapd m5, m0
	63	+ addpd m0, [sub_tab]
	64	+
	65	+ lea off2q, [lenq*4 - mmsize/2]
	66	+ sub lenq, mmsize/4 ; avoid overwriting
	67	+ xor off1q, off1q
	68	+
	69	+ cmp lenq, mmsize/4
	70	+ jl .scalar_o
	71	+
	72	+%if cpuflag(avx2)
	73	+ movapd m7, [dec_tab_avx2]
	74	+%else
	75	+ movapd m7, [dec_tab_sse2]
	76	+%endif
	77	+
	78	+.loop_o:
	79	+ movapd m1, m6
	80	+ mulpd m2, m0, m0
	81	+ subpd m1, m2
	82	+%if cpuflag(avx2)
	83	+ vpermpd m2, m1, q0123
	84	+%else
	85	+ shufpd m2, m1, m1, 01b
	86	+%endif
	87	+
	88	+ cvtdq2pd m3, [dataq + off1q]
	89	+ cvtdq2pd m4, [dataq + off2q]
	90	+
	91	+ mulpd m1, m3
	92	+ mulpd m2, m4
	93	+
	94	+ movupd [outq + off1q*2], m1
	95	+ movupd [outq + off2q*2], m2
	96	+
	97	+ addpd m0, m7
	98	+ add off1q, mmsize/2
	99	+ sub off2q, mmsize/2
	100	+ sub lenq, mmsize/4
	101	+ jg .loop_o
	102	+
	103	+ add lend, (mmsize/4 - 1)
	104	+ cmp lend, 0
	105	+ je .end
	106	+ sub lenq, (mmsize/4 - 1)
	107	+
	108	+.scalar_o:
	109	+ movapd xm7, [dec_tab_scalar]
	110	+ subpd xm0, xm7
	111	+
	112	+ ; Set offsets
	113	+ add off2q, (mmsize/4) + 4*cpuflag(avx2)
	114	+ add lenq, mmsize/4 - 2
	115	+
	116	+.loop_o_scalar:
	117	+ movapd xm1, xm6
	118	+ mulpd xm2, xm0, xm0
	119	+ subpd xm1, xm2
	120	+
	121	+ cvtdq2pd xm3, [dataq + off1q - (mmsize/4) + 4*cpuflag(avx2)]
	122	+ cvtdq2pd xm4, [dataq + off2q - (mmsize/4) + 4*cpuflag(avx2)]
	123	+
	124	+ mulpd xm3, xm1
	125	+ mulpd xm4, xm1
	126	+
	127	+ movhpd [outq + off1q*2], xm3
	128	+ movhpd [outq + off2q*2], xm4
	129	+
	130	+ addpd xm0, xm7
	131	+
	132	+ add off1q, 4
	133	+ sub off2q, 4
	134	+
	135	+ sub lenq, 2
	136	+ jg .loop_o_scalar
	137	+ RET
	138	+
	139	+.even:
	140	+%if cpuflag(avx2)
	141	+ addpd m0, [seq_tab_avx2]
	142	+%else
	143	+ addpd m0, [seq_tab_sse2]
	144	+%endif
	145	+
	146	+ mov off1d, lend
	147	+ shr off1d, 1
	148	+ movd xm1, off1d
	149	+ cvtdq2pd xm1, xm1 ; len/2
	150	+%if cpuflag(avx2)
	151	+ vbroadcastsd m1, xm1
	152	+%else
	153	+ shufpd m1, m1, 00b
	154	+%endif
	155	+ subpd m0, m1
	156	+
	157	+%if cpuflag(avx2)
	158	+ movapd m7, [add_tab_avx2]
	159	+%else
	160	+ movapd m7, [add_tab_sse2]
	161	+%endif
	162	+
	163	+ lea off2q, [lenq*2]
	164	+ lea off1q, [lenq*2 - mmsize/2]
	165	+ sub lenq, mmsize/4
	166	+
	167	+ cmp lenq, mmsize/4
	168	+ jl .scalar_e
	169	+
	170	+.loop_e:
	171	+ movapd m1, m6
	172	+ mulpd m2, m0, m0
	173	+ subpd m1, m2
	174	+%if cpuflag(avx2)
	175	+ vpermpd m2, m1, q0123
	176	+%else
	177	+ shufpd m2, m1, m1, 01b
	178	+%endif
	179	+
	180	+ cvtdq2pd m3, [dataq + off1q]
	181	+ cvtdq2pd m4, [dataq + off2q]
	182	+
	183	+ mulpd m1, m3
	184	+ mulpd m2, m4
	185	+
	186	+ movupd [outq + off1q*2], m1
	187	+ movupd [outq + off2q*2], m2
	188	+
	189	+ addpd m0, m7
	190	+ add off2q, mmsize/2
	191	+ sub off1q, mmsize/2
	192	+ sub lenq, mmsize/4
	193	+ jge .loop_e
	194	+
	195	+.scalar_e:
	196	+ subpd m0, m7
	197	+ movapd m7, [dec_tab_scalar]
	198	+ subpd m0, m7
	199	+ subpd m0, m7
	200	+ subpd m0, m7
	201	+
	202	+ add off1q, (mmsize/2)
	203	+ sub off2q, (mmsize/2) - 4 - 8*cpuflag(avx2)
	204	+
	205	+ addpd xm0, [sub_tab]
	206	+
	207	+.loop_e_scalar:
	208	+ movapd xm1, xm6
	209	+ mulpd xm2, xm0, xm0
	210	+ subpd xm1, xm2
	211	+
	212	+ cvtdq2pd m3, [dataq + off1q - 4]
	213	+ cvtdq2pd m4, [dataq + off2q - 4]
	214	+
	215	+ mulpd m3, m1
	216	+ mulpd m4, m1
	217	+
	218	+ movhpd [outq + off1q*2], xm3
	219	+ movhpd [outq + off2q*2], xm4
	220	+
	221	+ subpd xm0, xm7
	222	+
	223	+ add off2q, 4
	224	+ sub off1q, 4
	225	+ jge .loop_e_scalar
	226	+ RET
	227	+
	228	+.one:
	229	+ xorpd xm0, xm0
	230	+ movhpd [outq], xm0
	231	+.end:
	232	+ RET
	233	+%endmacro
	234	+
	235	+INIT_XMM sse2
	236	+APPLY_WELCH_FN
	237	+
	238	+%if HAVE_AVX2_EXTERNAL
	239	+INIT_YMM avx2
	240	+APPLY_WELCH_FN
	241	+%endif

libavcodec/x86/lpc_init.crenamed from libavcodec/x86/lpc.c+15−57

@@ -20,65 +20,19 @@ · 19 unmodified lines
20	20	*/
21	21
22	22	#include "libavutil/attributes.h"
23		−#include "libavutil/cpu.h"
24		−#include "libavutil/mem_internal.h"
25	23	#include "libavutil/x86/asm.h"
26	24	#include "libavutil/x86/cpu.h"
27	25	#include "libavcodec/lpc.h"
28	26
	27	+void ff_lpc_apply_welch_window_sse2(const int32_t *data, int len,
	28	+ double *w_data);
	29	+void ff_lpc_apply_welch_window_avx2(const int32_t *data, int len,
	30	+ double *w_data);
	31	+
29	32	DECLARE_ASM_CONST(16, double, pd_1)[2] = { 1.0, 1.0 };
30		−DECLARE_ASM_CONST(16, double, pd_2)[2] = { 2.0, 2.0 };
31	33
32	34	#if HAVE_SSE2_INLINE
33	35
34		−static void lpc_apply_welch_window_sse2(const int32_t *data, int len,
35		− double *w_data)
36		−{
37		− double c = 2.0 / (len-1.0);
38		− int n2 = len>>1;
39		− x86_reg i = -n2*sizeof(int32_t);
40		− x86_reg j = n2*sizeof(int32_t);
41		− __asm__ volatile(
42		− "movsd %4, %%xmm7 \n\t"
43		− "movapd "MANGLE(pd_1)", %%xmm6 \n\t"
44		− "movapd "MANGLE(pd_2)", %%xmm5 \n\t"
45		− "movlhps %%xmm7, %%xmm7 \n\t"
46		− "subpd %%xmm5, %%xmm7 \n\t"
47		− "addsd %%xmm6, %%xmm7 \n\t"
48		− "test $1, %5 \n\t"
49		− "jz 2f \n\t"
50		−#define WELCH(MOVPD, offset)\
51		− "1: \n\t"\
52		− "movapd %%xmm7, %%xmm1 \n\t"\
53		− "mulpd %%xmm1, %%xmm1 \n\t"\
54		− "movapd %%xmm6, %%xmm0 \n\t"\
55		− "subpd %%xmm1, %%xmm0 \n\t"\
56		− "pshufd $0x4e, %%xmm0, %%xmm1 \n\t"\
57		− "cvtpi2pd (%3,%0), %%xmm2 \n\t"\
58		− "cvtpi2pd "#offset"*4(%3,%1), %%xmm3 \n\t"\
59		− "mulpd %%xmm0, %%xmm2 \n\t"\
60		− "mulpd %%xmm1, %%xmm3 \n\t"\
61		− "movapd %%xmm2, (%2,%0,2) \n\t"\
62		− MOVPD" %%xmm3, "#offset"*8(%2,%1,2) \n\t"\
63		− "subpd %%xmm5, %%xmm7 \n\t"\
64		− "sub $8, %1 \n\t"\
65		− "add $8, %0 \n\t"\
66		− "jl 1b \n\t"\
67		−
68		− WELCH("movupd", -1)
69		− "jmp 3f \n\t"
70		− "2: \n\t"
71		− WELCH("movapd", -2)
72		− "3: \n\t"
73		− :"+&r"(i), "+&r"(j)
74		− :"r"(w_data+n2), "r"(data+n2), "m"(c), "r"(len)
75		− NAMED_CONSTRAINTS_ARRAY_ADD(pd_1,pd_2)
76		− XMM_CLOBBERS_ONLY("%xmm0", "%xmm1", "%xmm2", "%xmm3",
77		− "%xmm5", "%xmm6", "%xmm7")
78		− );
79		−#undef WELCH
80		−}
81		−
82	36	static void lpc_compute_autocorr_sse2(const double *data, int len, int lag,
83	37	double *autoc)
84	38	{
@@ -151,12 +105,16 @@ static void lpc_compute_autocorr_sse2(const double *data, int len, int lag, · 66 unmodified lines
151	105
152	106	av_cold void ff_lpc_init_x86(LPCContext *c)
153	107	{
154		−#if HAVE_SSE2_INLINE
155	108	int cpu_flags = av_get_cpu_flags();
156	109
157		− if (INLINE_SSE2_SLOW(cpu_flags)) {
158		− c->lpc_apply_welch_window = lpc_apply_welch_window_sse2;
159		− c->lpc_compute_autocorr = lpc_compute_autocorr_sse2;
160		− }
161		−#endif /* HAVE_SSE2_INLINE */
	110	+#if HAVE_SSE2_INLINE
	111	+ if (INLINE_SSE2_SLOW(cpu_flags))
	112	+ c->lpc_compute_autocorr = lpc_compute_autocorr_sse2;
	113	+#endif
	114	+
	115	+ if (EXTERNAL_SSE2(cpu_flags))
	116	+ c->lpc_apply_welch_window = ff_lpc_apply_welch_window_sse2;
	117	+
	118	+ if (EXTERNAL_AVX2(cpu_flags))
	119	+ c->lpc_apply_welch_window = ff_lpc_apply_welch_window_avx2;
162	120	}

tests/checkasm/Makefile+1−0

@@ -11,6 +11,7 @@ AVCODECOBJS-$(CONFIG_H264QPEL) += h264qpel.o · 10 unmodified lines
11	11	AVCODECOBJS-$(CONFIG_IDCTDSP) += idctdsp.o
12	12	AVCODECOBJS-$(CONFIG_LLVIDDSP) += llviddsp.o
13	13	AVCODECOBJS-$(CONFIG_LLVIDENCDSP) += llviddspenc.o
	14	+AVCODECOBJS-$(CONFIG_LPC) += lpc.o
14	15	AVCODECOBJS-$(CONFIG_ME_CMP) += motion.o
15	16	AVCODECOBJS-$(CONFIG_VC1DSP) += vc1dsp.o
16	17	AVCODECOBJS-$(CONFIG_VP8DSP) += vp8dsp.o

tests/checkasm/checkasm.c+3−0

@@ -135,6 +135,9 @@ static const struct { · 134 unmodified lines
135	135	#if CONFIG_LLVIDENCDSP
136	136	{ "llviddspenc", checkasm_check_llviddspenc },
137	137	#endif
	138	+ #if CONFIG_LPC
	139	+ { "lpc", checkasm_check_lpc },
	140	+ #endif
138	141	#if CONFIG_ME_CMP
139	142	{ "motion", checkasm_check_motion },
140	143	#endif

tests/checkasm/checkasm.h+1−0

@@ -68,6 +68,7 @@ void checkasm_check_idctdsp(void); · 67 unmodified lines
68	68	void checkasm_check_jpeg2000dsp(void);
69	69	void checkasm_check_llviddsp(void);
70	70	void checkasm_check_llviddspenc(void);
	71	+void checkasm_check_lpc(void);
71	72	void checkasm_check_motion(void);
72	73	void checkasm_check_nlmeans(void);
73	74	void checkasm_check_opusdsp(void);

tests/checkasm/lpc.c+80−0

@@ -0,0 +1,80 @@
	1	+/*
	2	+ * This file is part of FFmpeg.
	3	+ *
	4	+ * FFmpeg is free software; you can redistribute it and/or modify
	5	+ * it under the terms of the GNU General Public License as published by
	6	+ * the Free Software Foundation; either version 2 of the License, or
	7	+ * (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
	12	+ * GNU General Public License for more details.
	13	+ *
	14	+ * You should have received a copy of the GNU General Public License along
	15	+ * with FFmpeg; if not, write to the Free Software Foundation, Inc.,
	16	+ * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
	17	+ */
	18	+
	19	+#include "libavutil/mem_internal.h"
	20	+
	21	+#include "libavcodec/lpc.h"
	22	+
	23	+#include "checkasm.h"
	24	+
	25	+#define randomize_int32(buf, len) \
	26	+ do { \
	27	+ for (int i = 0; i < len; i++) { \
	28	+ int32_t f = rnd() >> 8; \
	29	+ buf[i] = f; \
	30	+ } \
	31	+ } while (0)
	32	+
	33	+#define EPS 0.005
	34	+
	35	+static void test_window(int len)
	36	+{
	37	+ LOCAL_ALIGNED(16, int32_t, src, [5000]);
	38	+ LOCAL_ALIGNED(16, double, dst0, [5000]);
	39	+ LOCAL_ALIGNED(16, double, dst1, [5000]);
	40	+
	41	+ declare_func(void, int32_t in, int len, double out);
	42	+
	43	+ randomize_int32(src, len);
	44	+
	45	+ call_ref(src, len, dst0);
	46	+ call_new(src, len, dst1);
	47	+
	48	+ if (!double_near_abs_eps_array(dst0, dst1, EPS, len))
	49	+ fail();
	50	+
	51	+ bench_new(src, len, dst1);
	52	+}
	53	+
	54	+void checkasm_check_lpc(void)
	55	+{
	56	+ LPCContext ctx;
	57	+ ff_lpc_init(&ctx, 32, 16, FF_LPC_TYPE_DEFAULT);
	58	+
	59	+ if (check_func(ctx.lpc_apply_welch_window, "apply_welch_window_even")) {
	60	+ for (int i = 0; i < 64; i += 2)
	61	+ test_window(i);
	62	+ }
	63	+ report("apply_welch_window_even");
	64	+
	65	+ if (check_func(ctx.lpc_apply_welch_window, "apply_welch_window_odd")) {
	66	+ for (int i = 1; i < 64; i += 2)
	67	+ test_window(i);
	68	+ }
	69	+ report("apply_welch_window_odd");
	70	+
	71	+ if (check_func(ctx.lpc_apply_welch_window, "apply_welch_window_4096"))
	72	+ test_window(4096);
	73	+ report("apply_welch_window_4096");
	74	+
	75	+ if (check_func(ctx.lpc_apply_welch_window, "apply_welch_window_4097"))
	76	+ test_window(4097);
	77	+ report("apply_welch_window_4097");
	78	+
	79	+ ff_lpc_end(&ctx);
	80	+}

x86/lpc: implement a new Welch windowing function

Keyboard shortcuts

General

Issues and pull requests

Issues, pull requests and Inbox

Repository code