diff options
Diffstat (limited to 'Transceiver52M/arch')
28 files changed, 2321 insertions, 0 deletions
diff --git a/Transceiver52M/arch/Makefile.am b/Transceiver52M/arch/Makefile.am new file mode 100644 index 0000000..14e6c82 --- /dev/null +++ b/Transceiver52M/arch/Makefile.am @@ -0,0 +1,8 @@ +include $(top_srcdir)/Makefile.common + +SUBDIRS = common +if ARCH_ARM +SUBDIRS += arm +else +SUBDIRS += x86 +endif diff --git a/Transceiver52M/arch/arm/Makefile.am b/Transceiver52M/arch/arm/Makefile.am new file mode 100644 index 0000000..89ffb32 --- /dev/null +++ b/Transceiver52M/arch/arm/Makefile.am @@ -0,0 +1,22 @@ +if ARCH_ARM_A15 +ARCH_FLAGS = -mfpu=neon-vfpv4 +else +ARCH_FLAGS = -mfpu=neon +endif + +AM_CFLAGS = -Wall $(ARCH_FLAGS) -std=gnu99 -I../common +AM_CCASFLAGS = $(ARCH_FLAGS) + +noinst_LTLIBRARIES = libarch.la + +libarch_la_LIBADD = $(top_builddir)/Transceiver52M/arch/common/libarch_common.la + +libarch_la_SOURCES = \ + convert.c \ + convert_neon.S \ + convolve.c \ + convolve_neon.S \ + scale.c \ + scale_neon.S \ + mult.c \ + mult_neon.S diff --git a/Transceiver52M/arch/arm/convert.c b/Transceiver52M/arch/arm/convert.c new file mode 100644 index 0000000..c94a3d7 --- /dev/null +++ b/Transceiver52M/arch/arm/convert.c @@ -0,0 +1,85 @@ +/* + * NEON type conversions + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include "convert.h" + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +void neon_convert_ps_si16_4n(short *, const float *, const float *, int); +void neon_convert_si16_ps_4n(float *, const short *, int); + +void convert_init(void) { +} + +/* 4*N 16-bit signed integer conversion with remainder */ +static void neon_convert_si16_ps(float *out, + const short *in, + int len) +{ + int start = len / 4 * 4; + + neon_convert_si16_ps_4n(out, in, len >> 2); + + for (int i = 0; i < len % 4; i++) + out[start + i] = (float) in[start + i]; +} + +/* 4*N 16-bit signed integer conversion with remainder */ +static void neon_convert_ps_si16(short *out, + const float *in, + const float *scale, + int len) +{ + int start = len / 4 * 4; + + neon_convert_ps_si16_4n(out, in, scale, len >> 2); + + for (int i = 0; i < len % 4; i++) + out[start + i] = (short) (in[start + i] * (*scale)); +} + +void convert_float_short(short *out, const float *in, float scale, int len) +{ +#ifdef HAVE_NEON + float q[4] = { scale, scale, scale, scale }; + + if (len % 4) + neon_convert_ps_si16(out, in, q, len); + else + neon_convert_ps_si16_4n(out, in, q, len >> 2); +#else + base_convert_float_short(out, in, scale, len); +#endif +} + +void convert_short_float(float *out, const short *in, int len) +{ +#ifdef HAVE_NEON + if (len % 4) + neon_convert_si16_ps(out, in, len); + else + neon_convert_si16_ps_4n(out, in, len >> 2); +#else + base_convert_short_float(out, in, len); +#endif +} diff --git a/Transceiver52M/arch/arm/convert_neon.S b/Transceiver52M/arch/arm/convert_neon.S new file mode 100644 index 0000000..842ed9f --- /dev/null +++ b/Transceiver52M/arch/arm/convert_neon.S @@ -0,0 +1,51 @@ +/* + * NEON type conversions + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + + .syntax unified + .text + .align 2 + .global neon_convert_ps_si16_4n + .type neon_convert_ps_si16_4n, %function +neon_convert_ps_si16_4n: + vld1.32 {q1}, [r2] +.loop_fltint: + vld1.64 {d0-d1}, [r1]! + vmul.f32 q0, q1 + vcvt.s32.f32 q2, q0 + vqmovn.s32 d0, q2 + vst1.64 {d0}, [r0]! + subs r3, #1 + bne .loop_fltint + bx lr + .size neon_convert_ps_si16_4n, .-neon_convert_ps_si16_4n + .text + .align 2 + .global neon_convert_si16_ps_4n + .type neon_convert_si16_ps_4n, %function +neon_convert_si16_ps_4n: +.loop_intflt: + vld1.64 {d0}, [r1]! + vmovl.s16 q1, d0 + vcvt.f32.s32 q0, q1 + vst1.64 {q0}, [r0]! + subs r2, #1 + bne .loop_intflt + bx lr + .size neon_convert_si16_ps_4n, .-neon_convert_si16_ps_4n + .section .note.GNU-stack,"",%progbits diff --git a/Transceiver52M/arch/arm/convolve.c b/Transceiver52M/arch/arm/convolve.c new file mode 100644 index 0000000..912d0c2 --- /dev/null +++ b/Transceiver52M/arch/arm/convolve.c @@ -0,0 +1,146 @@ +/* + * NEON Convolution + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include <stdio.h> + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +/* Forward declarations from base implementation */ +int _base_convolve_real(float *x, int x_len, + float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +int _base_convolve_complex(float *x, int x_len, + float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +int bounds_check(int x_len, int h_len, int y_len, + int start, int len, int step); + +#ifdef HAVE_NEON +/* Calls into NEON assembler */ +void neon_conv_real4(float *x, float *h, float *y, int len); +void neon_conv_real8(float *x, float *h, float *y, int len); +void neon_conv_real12(float *x, float *h, float *y, int len); +void neon_conv_real16(float *x, float *h, float *y, int len); +void neon_conv_real20(float *x, float *h, float *y, int len); +void mac_cx_neon4(float *x, float *h, float *y, int len); + +/* Complex-complex convolution */ +static void neon_conv_cmplx_4n(float *x, float *h, float *y, int h_len, int len) +{ + for (int i = 0; i < len; i++) + mac_cx_neon4(&x[2 * i], h, &y[2 * i], h_len >> 2); +} +#endif + +/* API: Initalize convolve module */ +void convolve_init(void) +{ + /* Stub */ + return; +} + +/* API: Aligned complex-real */ +int convolve_real(float *x, int x_len, + float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset) +{ + void (*conv_func)(float *, float *, float *, int) = NULL; + + if (bounds_check(x_len, h_len, y_len, start, len, step) < 0) + return -1; + + memset(y, 0, len * 2 * sizeof(float)); + +#ifdef HAVE_NEON + if (step <= 4) { + switch (h_len) { + case 4: + conv_func = neon_conv_real4; + break; + case 8: + conv_func = neon_conv_real8; + break; + case 12: + conv_func = neon_conv_real12; + break; + case 16: + conv_func = neon_conv_real16; + break; + case 20: + conv_func = neon_conv_real20; + break; + } + } +#endif + if (conv_func) { + conv_func(&x[2 * (-(h_len - 1) + start)], + h, y, len); + } else { + _base_convolve_real(x, x_len, + h, h_len, + y, y_len, + start, len, step, offset); + } + + return len; +} + + +/* API: Aligned complex-complex */ +int convolve_complex(float *x, int x_len, + float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset) +{ + void (*conv_func)(float *, float *, float *, int, int) = NULL; + + if (bounds_check(x_len, h_len, y_len, start, len, step) < 0) + return -1; + + memset(y, 0, len * 2 * sizeof(float)); + +#ifdef HAVE_NEON + if (step <= 4 && !(h_len % 4)) + conv_func = neon_conv_cmplx_4n; +#endif + if (conv_func) { + conv_func(&x[2 * (-(h_len - 1) + start)], + h, y, h_len, len); + } else { + _base_convolve_complex(x, x_len, + h, h_len, + y, y_len, + start, len, step, offset); + } + + return len; +} diff --git a/Transceiver52M/arch/arm/convolve_neon.S b/Transceiver52M/arch/arm/convolve_neon.S new file mode 100644 index 0000000..637d150 --- /dev/null +++ b/Transceiver52M/arch/arm/convolve_neon.S @@ -0,0 +1,277 @@ +/* + * NEON Convolution + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + + .syntax unified + .text + .align 2 + .global neon_conv_real4 + .type neon_conv_real4, %function +neon_conv_real4: + push {r4, lr} + vpush {q4-q7} + vld2.32 {q0-q1}, [r1] + ldr r4, =8 +.neon_conv_loop4: + vld2.32 {q2-q3}, [r0], r4 + vmul.f32 q4, q2, q0 + vmul.f32 q5, q3, q0 + vpadd.f32 d12, d8, d9 + vpadd.f32 d13, d10, d11 + vpadd.f32 d14, d12, d13 + vst1.64 {d14}, [r2]! + subs r3, r3, #1 + bne .neon_conv_loop4 + vpop {q4-q7} + pop {r4, pc} + .size neon_conv_real4, .-neon_conv_real4 + .align 2 + .p2align 4,,15 + .global neon_conv_real8 + .type neon_conv_real8, %function +neon_conv_real8: + push {r4-r5, lr} + vpush {q4-q7} + vld2.32 {q0-q1}, [r1]! + vld2.32 {q2-q3}, [r1] + add r4, r0, #32 + ldr r5, =8 +.neon_conv_loop8: + vld2.32 {q4-q5}, [r0], r5 + vld2.32 {q6-q7}, [r4], r5 + vmul.f32 q8, q4, q0 + vmul.f32 q9, q5, q0 + vmul.f32 q10, q6, q2 + vmul.f32 q11, q7, q2 + + vadd.f32 q12, q8, q10 + vadd.f32 q13, q9, q11 + + vpadd.f32 d22, d24, d25 + vpadd.f32 d23, d26, d27 + vpadd.f32 d24, d22, d23 + vst1.64 {d24}, [r2]! + subs r3, r3, #1 + bne .neon_conv_loop8 + vpop {q4-q7} + pop {r4-r5, pc} + .size neon_conv_real8, .-neon_conv_real8 + .align 2 + .global neon_conv_real12 + .type neon_conv_real12, %function +neon_conv_real12: + push {r4-r6, lr} + vpush {q4-q7} + vld2.32 {q0-q1}, [r1]! + vld2.32 {q2-q3}, [r1]! + vld2.32 {q4-q5}, [r1]! + add r4, r0, #32 + add r5, r0, #64 + ldr r6, =8 +.neon_conv_loop12: + vld2.32 {q6-q7}, [r0], r6 + vld2.32 {q8-q9}, [r4], r6 + vld2.32 {q10-q11}, [r5], r6 +#ifdef HAVE_NEON_FMA + vfma.f32 q1, q6, q0 + vfma.f32 q3, q7, q0 + vfma.f32 q1, q8, q2 + vfma.f32 q3, q9, q2 + vfma.f32 q1, q10, q4 + vfma.f32 q3, q11, q4 +#else + vmul.f32 q12, q6, q0 + vmul.f32 q13, q7, q0 + vmul.f32 q14, q8, q2 + vmul.f32 q15, q9, q2 + vmul.f32 q1, q10, q4 + vmul.f32 q3, q11, q4 + + vadd.f32 q5, q12, q14 + vadd.f32 q6, q13, q15 + vadd.f32 q1, q5, q1 + vadd.f32 q3, q6, q3 +#endif + vpadd.f32 d2, d2, d3 + vpadd.f32 d3, d6, d7 + vpadd.f32 d6, d2, d3 + vst1.64 {d6}, [r2]! + subs r3, r3, #1 + bne .neon_conv_loop12 + vpop {q4-q7} + pop {r4-r6, pc} + .size neon_conv_real12, .-neon_conv_real12 + .align 2 + .global neon_conv_real16 + .type neon_conv_real16, %function +neon_conv_real16: + push {r4-r7, lr} + vpush {q4-q7} + vld2.32 {q0-q1}, [r1]! + vld2.32 {q2-q3}, [r1]! + vld2.32 {q4-q5}, [r1]! + vld2.32 {q6-q7}, [r1] + add r4, r0, #32 + add r5, r0, #64 + add r6, r0, #96 + ldr r7, =8 +.neon_conv_loop16: + vld2.32 {q8-q9}, [r0], r7 + vld2.32 {q10-q11}, [r4], r7 + vld2.32 {q12-q13}, [r5], r7 + vld2.32 {q14-q15}, [r6], r7 +#ifdef HAVE_NEON_FMA + vmul.f32 q1, q8, q0 + vmul.f32 q3, q9, q0 + vfma.f32 q1, q10, q2 + vfma.f32 q3, q11, q2 + vfma.f32 q1, q12, q4 + vfma.f32 q3, q13, q4 + vfma.f32 q1, q14, q6 + vfma.f32 q3, q15, q6 +#else + vmul.f32 q1, q8, q0 + vmul.f32 q3, q9, q0 + vmul.f32 q5, q10, q2 + vmul.f32 q7, q11, q2 + vmul.f32 q8, q12, q4 + vmul.f32 q9, q13, q4 + vmul.f32 q10, q14, q6 + vmul.f32 q11, q15, q6 + + vadd.f32 q1, q1, q5 + vadd.f32 q3, q3, q7 + vadd.f32 q5, q8, q10 + vadd.f32 q7, q9, q11 + vadd.f32 q1, q1, q5 + vadd.f32 q3, q3, q7 +#endif + vpadd.f32 d2, d2, d3 + vpadd.f32 d3, d6, d7 + vpadd.f32 d6, d2, d3 + vst1.64 {d6}, [r2]! + subs r3, r3, #1 + bne .neon_conv_loop16 + vpop {q4-q7} + pop {r4-r7, pc} + .size neon_conv_real16, .-neon_conv_real16 + .align 2 + .global neon_conv_real20 + .type neon_conv_real20, %function +neon_conv_real20: + push {r4-r8, lr} + vpush {q4-q7} + vld2.32 {q0-q1}, [r1]! + vld2.32 {q2-q3}, [r1]! + vld2.32 {q4-q5}, [r1]! + vld2.32 {q6-q7}, [r1]! + vld2.32 {q8-q9}, [r1] + add r4, r0, #32 + add r5, r0, #64 + add r6, r0, #96 + add r7, r0, #128 + ldr r8, =8 +.neon_conv_loop20: + vld2.32 {q10-q11}, [r0], r8 + vld2.32 {q12-q13}, [r4], r8 + vld2.32 {q14-q15}, [r5], r8 +#ifdef HAVE_NEON_FMA + vmul.f32 q1, q10, q0 + vfma.f32 q1, q12, q2 + vfma.f32 q1, q14, q4 + vmul.f32 q3, q11, q0 + vfma.f32 q3, q13, q2 + vfma.f32 q3, q15, q4 + + vld2.32 {q12-q13}, [r6], r8 + vld2.32 {q14-q15}, [r7], r8 + + vfma.f32 q1, q12, q6 + vfma.f32 q3, q13, q6 + vfma.f32 q1, q14, q8 + vfma.f32 q3, q15, q8 +#else + vmul.f32 q1, q10, q0 + vmul.f32 q3, q12, q2 + vmul.f32 q5, q14, q4 + vmul.f32 q7, q11, q0 + vmul.f32 q9, q13, q2 + vmul.f32 q10, q15, q4 + vadd.f32 q1, q1, q3 + vadd.f32 q3, q7, q9 + vadd.f32 q9, q1, q5 + vadd.f32 q10, q3, q10 + + vld2.32 {q12-q13}, [r6], r8 + vld2.32 {q14-q15}, [r7], r8 + + vmul.f32 q1, q12, q6 + vmul.f32 q3, q13, q6 + vmul.f32 q5, q14, q8 + vmul.f32 q7, q15, q8 + vadd.f32 q12, q1, q9 + vadd.f32 q14, q3, q10 + vadd.f32 q1, q12, q5 + vadd.f32 q3, q14, q7 +#endif + vpadd.f32 d2, d2, d3 + vpadd.f32 d3, d6, d7 + vpadd.f32 d6, d2, d3 + vst1.64 {d6}, [r2]! + subs r3, r3, #1 + bne .neon_conv_loop20 + vpop {q4-q7} + pop {r4-r8, pc} + .size neon_conv_real20, .-neon_conv_real20 + .align 2 + .global mac_cx_neon4 + .type mac_cx_neon4, %function +mac_cx_neon4: + push {r4, lr} + ldr r4, =32 + veor q14, q14 + veor q15, q15 +.neon_conv_loop_mac4: + vld2.32 {q0-q1}, [r0], r4 + vld2.32 {q2-q3}, [r1]! + + vmul.f32 q10, q0, q2 + vmul.f32 q11, q1, q3 + vmul.f32 q12, q0, q3 + vmul.f32 q13, q2, q1 + vsub.f32 q8, q10, q11 + vadd.f32 q9, q12, q13 + + vadd.f32 q14, q8 + vadd.f32 q15, q9 + subs r3, #1 + bne .neon_conv_loop_mac4 + + vld1.64 d0, [r2] + vpadd.f32 d28, d28, d29 + vpadd.f32 d30, d30, d31 + vpadd.f32 d1, d28, d30 + vadd.f32 d1, d0 + vst1.64 d1, [r2] + pop {r4, pc} + .size mac_cx_neon4, .-mac_cx_neon4 + .section .note.GNU-stack,"",%progbits diff --git a/Transceiver52M/arch/arm/mult.c b/Transceiver52M/arch/arm/mult.c new file mode 100644 index 0000000..245be50 --- /dev/null +++ b/Transceiver52M/arch/arm/mult.c @@ -0,0 +1,56 @@ +/* + * NEON scaling + * Copyright (C) 2012,2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include <mult.h> + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +void neon_cmplx_mul_4n(float *, float *, float *, int); + +static void cmplx_mul_ps(float *out, float *a, float *b, int len) +{ + float ai, aq, bi, bq; + + for (int i = 0; i < len; i++) { + ai = a[2 * i + 0]; + aq = a[2 * i + 1]; + + bi = b[2 * i + 0]; + bq = b[2 * i + 1]; + + out[2 * i + 0] = ai * bi - aq * bq; + out[2 * i + 1] = ai * bq + aq * bi; + } +} + +void mul_complex(float *out, float *a, float *b, int len) +{ +#ifdef HAVE_NEON + if (len % 4) + cmplx_mul_ps(out, a, b, len); + else + neon_cmplx_mul_4n(out, a, b, len >> 2); +#else + cmplx_mul_ps(out, a, b, len); +#endif +} diff --git a/Transceiver52M/arch/arm/mult_neon.S b/Transceiver52M/arch/arm/mult_neon.S new file mode 100644 index 0000000..162846e --- /dev/null +++ b/Transceiver52M/arch/arm/mult_neon.S @@ -0,0 +1,42 @@ +/* + * NEON complex multiplication + * Copyright (C) 2012,2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + + .syntax unified + .text + .align 2 + .global neon_cmplx_mul_4n + .type neon_cmplx_mul_4n, %function +neon_cmplx_mul_4n: + vpush {q4-q7} +.loop_mul: + vld2.32 {q0-q1}, [r1]! + vld2.32 {q2-q3}, [r2]! + vmul.f32 q4, q0, q2 + vmul.f32 q5, q1, q3 + vmul.f32 q6, q0, q3 + vmul.f32 q7, q2, q1 + vsub.f32 q8, q4, q5 + vadd.f32 q9, q6, q7 + vst2.32 {q8-q9}, [r0]! + subs r3, #1 + bne .loop_mul + vpop {q4-q7} + bx lr + .size neon_cmplx_mul_4n, .-neon_cmplx_mul_4n + .section .note.GNU-stack,"",%progbits diff --git a/Transceiver52M/arch/arm/scale.c b/Transceiver52M/arch/arm/scale.c new file mode 100644 index 0000000..2de13ff --- /dev/null +++ b/Transceiver52M/arch/arm/scale.c @@ -0,0 +1,56 @@ +/* + * NEON scaling + * Copyright (C) 2012,2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include <scale.h> + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +void neon_scale_4n(float *, float *, float *, int); + +static void scale_ps(float *out, float *in, float *scale, int len) +{ + float ai, aq, bi, bq; + + bi = scale[0]; + bq = scale[1]; + + for (int i = 0; i < len; i++) { + ai = in[2 * i + 0]; + aq = in[2 * i + 1]; + + out[2 * i + 0] = ai * bi - aq * bq; + out[2 * i + 1] = ai * bq + aq * bi; + } +} + +void scale_complex(float *out, float *in, float* scale, int len) +{ +#ifdef HAVE_NEON + if (len % 4) + scale_ps(out, in, scale, len); + else + neon_scale_4n(in, scale, out, len >> 2); +#else + scale_ps(out, in, scale, len); +#endif +} diff --git a/Transceiver52M/arch/arm/scale_neon.S b/Transceiver52M/arch/arm/scale_neon.S new file mode 100644 index 0000000..a66fbe5 --- /dev/null +++ b/Transceiver52M/arch/arm/scale_neon.S @@ -0,0 +1,50 @@ +/* + * ARM NEON Scaling + * Copyright (C) 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + + .syntax unified + .text + .align 2 + .global neon_scale_4n + .type neon_scale_4n, %function +neon_scale_4n: + push {r4, lr} + ldr r4, =32 + + vld1.64 d0, [r1] + vmov.32 s4, s1 + vmov.32 s1, s0 + vmov.64 d1, d0 + vmov.32 s5, s4 + vmov.64 d3, d2 +.loop_mul_const: + vld2.32 {q2-q3}, [r0], r4 + + vmul.f32 q8, q0, q2 + vmul.f32 q9, q1, q3 + vmul.f32 q10, q0, q3 + vmul.f32 q11, q1, q2 + vsub.f32 q8, q8, q9 + vadd.f32 q9, q10, q11 + + vst2.32 {q8-q9}, [r2]! + subs r3, #1 + bne .loop_mul_const + pop {r4, pc} + .size neon_scale_4n, .-neon_scale_4n + .section .note.GNU-stack,"",%progbits diff --git a/Transceiver52M/arch/common/Makefile.am b/Transceiver52M/arch/common/Makefile.am new file mode 100644 index 0000000..6b37906 --- /dev/null +++ b/Transceiver52M/arch/common/Makefile.am @@ -0,0 +1,15 @@ +AM_CFLAGS = -Wall -std=gnu99 + +noinst_LTLIBRARIES = libarch_common.la + +noinst_HEADERS = \ + convolve.h \ + convert.h \ + scale.h \ + mult.h \ + fft.h + +libarch_common_la_SOURCES = \ + convolve_base.c \ + convert_base.c \ + fft.c diff --git a/Transceiver52M/arch/common/convert.h b/Transceiver52M/arch/common/convert.h new file mode 100644 index 0000000..73402b0 --- /dev/null +++ b/Transceiver52M/arch/common/convert.h @@ -0,0 +1,15 @@ +#ifndef _CONVERT_H_ +#define _CONVERT_H_ + +void convert_float_short(short *out, const float *in, float scale, int len); + +void convert_short_float(float *out, const short *in, int len); + +void base_convert_float_short(short *out, const float *in, + float scale, int len); + +void base_convert_short_float(float *out, const short *in, int len); + +void convert_init(void); + +#endif /* _CONVERT_H_ */ diff --git a/Transceiver52M/arch/common/convert_base.c b/Transceiver52M/arch/common/convert_base.c new file mode 100644 index 0000000..5251fb8 --- /dev/null +++ b/Transceiver52M/arch/common/convert_base.c @@ -0,0 +1,34 @@ +/* + * Conversion + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include "convert.h" + +void base_convert_float_short(short *out, const float *in, + float scale, int len) +{ + for (int i = 0; i < len; i++) + out[i] = in[i] * scale; +} + +void base_convert_short_float(float *out, const short *in, int len) +{ + for (int i = 0; i < len; i++) + out[i] = in[i]; +} + diff --git a/Transceiver52M/arch/common/convolve.h b/Transceiver52M/arch/common/convolve.h new file mode 100644 index 0000000..43db577 --- /dev/null +++ b/Transceiver52M/arch/common/convolve.h @@ -0,0 +1,32 @@ +#ifndef _CONVOLVE_H_ +#define _CONVOLVE_H_ + +void *convolve_h_alloc(int num); + +int convolve_real(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +int convolve_complex(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +int base_convolve_real(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +int base_convolve_complex(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +void convolve_init(void); + +#endif /* _CONVOLVE_H_ */ diff --git a/Transceiver52M/arch/common/convolve_base.c b/Transceiver52M/arch/common/convolve_base.c new file mode 100644 index 0000000..71453a1 --- /dev/null +++ b/Transceiver52M/arch/common/convolve_base.c @@ -0,0 +1,156 @@ +/* + * Convolution + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include <stdio.h> + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +/* Base multiply and accumulate complex-real */ +static void mac_real(const float *x, const float *h, float *y) +{ + y[0] += x[0] * h[0]; + y[1] += x[1] * h[0]; +} + +/* Base multiply and accumulate complex-complex */ +static void mac_cmplx(const float *x, const float *h, float *y) +{ + y[0] += x[0] * h[0] - x[1] * h[1]; + y[1] += x[0] * h[1] + x[1] * h[0]; +} + +/* Base vector complex-complex multiply and accumulate */ +static void mac_real_vec_n(const float *x, const float *h, float *y, + int len, int step, int offset) +{ + for (int i = offset; i < len; i += step) + mac_real(&x[2 * i], &h[2 * i], y); +} + +/* Base vector complex-complex multiply and accumulate */ +static void mac_cmplx_vec_n(const float *x, const float *h, float *y, + int len, int step, int offset) +{ + for (int i = offset; i < len; i += step) + mac_cmplx(&x[2 * i], &h[2 * i], y); +} + +/* Base complex-real convolution */ +int _base_convolve_real(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset) +{ + for (int i = 0; i < len; i++) { + mac_real_vec_n(&x[2 * (i - (h_len - 1) + start)], + h, + &y[2 * i], h_len, + step, offset); + } + + return len; +} + +/* Base complex-complex convolution */ +int _base_convolve_complex(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset) +{ + for (int i = 0; i < len; i++) { + mac_cmplx_vec_n(&x[2 * (i - (h_len - 1) + start)], + h, + &y[2 * i], + h_len, step, offset); + } + + return len; +} + +/* Buffer validity checks */ +int bounds_check(int x_len, int h_len, int y_len, + int start, int len, int step) +{ + if ((x_len < 1) || (h_len < 1) || + (y_len < 1) || (len < 1) || (step < 1)) { + fprintf(stderr, "Convolve: Invalid input\n"); + return -1; + } + + if ((start + len > x_len) || (len > y_len) || (x_len < h_len)) { + fprintf(stderr, "Convolve: Boundary exception\n"); + fprintf(stderr, "start: %i, len: %i, x: %i, h: %i, y: %i\n", + start, len, x_len, h_len, y_len); + return -1; + } + + return 0; +} + +/* API: Non-aligned (no SSE) complex-real */ +int base_convolve_real(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset) +{ + if (bounds_check(x_len, h_len, y_len, start, len, step) < 0) + return -1; + + memset(y, 0, len * 2 * sizeof(float)); + + return _base_convolve_real(x, x_len, + h, h_len, + y, y_len, + start, len, step, offset); +} + +/* API: Non-aligned (no SSE) complex-complex */ +int base_convolve_complex(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset) +{ + if (bounds_check(x_len, h_len, y_len, start, len, step) < 0) + return -1; + + memset(y, 0, len * 2 * sizeof(float)); + + return _base_convolve_complex(x, x_len, + h, h_len, + y, y_len, + start, len, step, offset); +} + +/* Aligned filter tap allocation */ +void *convolve_h_alloc(int len) +{ +#ifdef HAVE_SSE3 + return memalign(16, len * 2 * sizeof(float)); +#else + return malloc(len * 2 * sizeof(float)); +#endif +} diff --git a/Transceiver52M/arch/common/fft.c b/Transceiver52M/arch/common/fft.c new file mode 100644 index 0000000..18b2de7 --- /dev/null +++ b/Transceiver52M/arch/common/fft.c @@ -0,0 +1,112 @@ +/* + * Fast Fourier transform + * + * Copyright (C) 2012 Tom Tsou <tom@tsou.cc> + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program 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 Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program; if not, see <http://www.gnu.org/licenses/>. + * See the COPYING file in the main directory for details. + */ + +#include <stdlib.h> +#include <string.h> +#include <assert.h> +#include <fftw3.h> + +#include "fft.h" + +struct fft_hdl { + float *fft_in; + float *fft_out; + int len; + fftwf_plan fft_plan; +}; + +/*! \brief Initialize FFT backend + * \param[in] reverse FFT direction + * \param[in] m FFT length + * \param[in] istride input stride count + * \param[in] ostride output stride count + * \param[in] in input buffer (FFTW aligned) + * \param[in] out output buffer (FFTW aligned) + * \param[in] ooffset initial offset into output buffer + * + * If the reverse is non-NULL, then an inverse FFT will be used. This is a + * wrapper for advanced non-contiguous FFTW usage. See FFTW documentation for + * further details. + * + * http://www.fftw.org/doc/Advanced-Complex-DFTs.html + * + * It is currently unknown how the offset of the output buffer affects FFTW + * memory alignment. + */ +struct fft_hdl *init_fft(int reverse, int m, int istride, int ostride, + float *in, float *out, int ooffset) +{ + int rank = 1; + int n[] = { m }; + int howmany = istride; + int idist = 1; + int odist = 1; + int *inembed = n; + int *onembed = n; + fftwf_complex *obuffer, *ibuffer; + + struct fft_hdl *hdl = (struct fft_hdl *) malloc(sizeof(struct fft_hdl)); + if (!hdl) + return NULL; + + int direction = FFTW_FORWARD; + if (reverse) + direction = FFTW_BACKWARD; + + ibuffer = (fftwf_complex *) in; + obuffer = (fftwf_complex *) out + ooffset; + + hdl->fft_in = in; + hdl->fft_out = out; + hdl->fft_plan = fftwf_plan_many_dft(rank, n, howmany, + ibuffer, inembed, istride, idist, + obuffer, onembed, ostride, odist, + direction, FFTW_MEASURE); + return hdl; +} + +void *fft_malloc(size_t size) +{ + return fftwf_malloc(size); +} + +void fft_free(void *ptr) +{ + free(ptr); +} + +/*! \brief Free FFT backend resources + */ +void free_fft(struct fft_hdl *hdl) +{ + fftwf_destroy_plan(hdl->fft_plan); + free(hdl); +} + +/*! \brief Run multiple DFT operations with the initialized plan + * \param[in] hdl handle to an intitialized fft struct + * + * Input and output buffers are configured with init_fft(). + */ +int cxvec_fft(struct fft_hdl *hdl) +{ + fftwf_execute(hdl->fft_plan); + return 0; +} diff --git a/Transceiver52M/arch/common/fft.h b/Transceiver52M/arch/common/fft.h new file mode 100644 index 0000000..fb7bede --- /dev/null +++ b/Transceiver52M/arch/common/fft.h @@ -0,0 +1,13 @@ +#ifndef _FFT_H_ +#define _FFT_H_ + +struct fft_hdl; + +struct fft_hdl *init_fft(int reverse, int m, int istride, int ostride, + float *in, float *out, int ooffset); +void *fft_malloc(size_t size); +void fft_free(void *ptr); +void free_fft(struct fft_hdl *hdl); +int cxvec_fft(struct fft_hdl *hdl); + +#endif /* _FFT_H_ */ diff --git a/Transceiver52M/arch/common/mult.h b/Transceiver52M/arch/common/mult.h new file mode 100644 index 0000000..4d96efb --- /dev/null +++ b/Transceiver52M/arch/common/mult.h @@ -0,0 +1,6 @@ +#ifndef _MULT_H_ +#define _MULT_H_ + +void mul_complex(float *out, float *a, float *b, int len); + +#endif /* _MULT_H_ */ diff --git a/Transceiver52M/arch/common/scale.h b/Transceiver52M/arch/common/scale.h new file mode 100644 index 0000000..da867e7 --- /dev/null +++ b/Transceiver52M/arch/common/scale.h @@ -0,0 +1,6 @@ +#ifndef _SCALE_H_ +#define _SCALE_H_ + +void scale_complex(float *out, float *in, float *scale, int len); + +#endif /* _SCALE_H_ */ diff --git a/Transceiver52M/arch/x86/Makefile.am b/Transceiver52M/arch/x86/Makefile.am new file mode 100644 index 0000000..f39dde5 --- /dev/null +++ b/Transceiver52M/arch/x86/Makefile.am @@ -0,0 +1,28 @@ +AM_CFLAGS = -Wall -std=gnu99 -I${srcdir}/../common + +noinst_LTLIBRARIES = libarch.la +noinst_LTLIBRARIES += libarch_sse_3.la +noinst_LTLIBRARIES += libarch_sse_4_1.la + +libarch_la_LIBADD = $(top_builddir)/Transceiver52M/arch/common/libarch_common.la + +# SSE 3 specific code +if HAVE_SSE3 +libarch_sse_3_la_SOURCES = \ + convert_sse_3.c \ + convolve_sse_3.c +libarch_sse_3_la_CFLAGS = $(AM_CFLAGS) -msse3 +libarch_la_LIBADD += libarch_sse_3.la +endif + +# SSE 4.1 specific code +if HAVE_SSE4_1 +libarch_sse_4_1_la_SOURCES = \ + convert_sse_4_1.c +libarch_sse_4_1_la_CFLAGS = $(AM_CFLAGS) -msse4.1 +libarch_la_LIBADD += libarch_sse_4_1.la +endif + +libarch_la_SOURCES = \ + convert.c \ + convolve.c diff --git a/Transceiver52M/arch/x86/convert.c b/Transceiver52M/arch/x86/convert.c new file mode 100644 index 0000000..07cdf59 --- /dev/null +++ b/Transceiver52M/arch/x86/convert.c @@ -0,0 +1,83 @@ +/* + * SSE type conversions + * Copyright (C) 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include "convert.h" +#include "convert_sse_3.h" +#include "convert_sse_4_1.h" + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +/* Architecture dependant function pointers */ +struct convert_cpu_context { + void (*convert_si16_ps_16n) (float *, const short *, int); + void (*convert_si16_ps) (float *, const short *, int); + void (*convert_scale_ps_si16_16n)(short *, const float *, float, int); + void (*convert_scale_ps_si16_8n)(short *, const float *, float, int); + void (*convert_scale_ps_si16)(short *, const float *, float, int); +}; + +static struct convert_cpu_context c; + +void convert_init(void) +{ + c.convert_scale_ps_si16_16n = base_convert_float_short; + c.convert_scale_ps_si16_8n = base_convert_float_short; + c.convert_scale_ps_si16 = base_convert_float_short; + c.convert_si16_ps_16n = base_convert_short_float; + c.convert_si16_ps = base_convert_short_float; + +#ifdef HAVE___BUILTIN_CPU_SUPPORTS +#ifdef HAVE_SSE4_1 + if (__builtin_cpu_supports("sse4.1")) { + c.convert_si16_ps_16n = &_sse_convert_si16_ps_16n; + c.convert_si16_ps = &_sse_convert_si16_ps; + } +#endif + +#ifdef HAVE_SSE3 + if (__builtin_cpu_supports("sse3")) { + c.convert_scale_ps_si16_16n = _sse_convert_scale_ps_si16_16n; + c.convert_scale_ps_si16_8n = _sse_convert_scale_ps_si16_8n; + c.convert_scale_ps_si16 = _sse_convert_scale_ps_si16; + } +#endif +#endif +} + +void convert_float_short(short *out, const float *in, float scale, int len) +{ + if (!(len % 16)) + c.convert_scale_ps_si16_16n(out, in, scale, len); + else if (!(len % 8)) + c.convert_scale_ps_si16_8n(out, in, scale, len); + else + c.convert_scale_ps_si16(out, in, scale, len); +} + +void convert_short_float(float *out, const short *in, int len) +{ + if (!(len % 16)) + c.convert_si16_ps_16n(out, in, len); + else + c.convert_si16_ps(out, in, len); +} diff --git a/Transceiver52M/arch/x86/convert_sse_3.c b/Transceiver52M/arch/x86/convert_sse_3.c new file mode 100644 index 0000000..255db67 --- /dev/null +++ b/Transceiver52M/arch/x86/convert_sse_3.c @@ -0,0 +1,107 @@ +/* + * SSE type conversions + * Copyright (C) 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include "convert_sse_3.h" + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_SSE3 +#include <xmmintrin.h> +#include <emmintrin.h> + +/* 8*N single precision floats scaled and converted to 16-bit signed integer */ +void _sse_convert_scale_ps_si16_8n(short *restrict out, + const float *restrict in, + float scale, int len) +{ + __m128 m0, m1, m2; + __m128i m4, m5; + + for (int i = 0; i < len / 8; i++) { + /* Load (unaligned) packed floats */ + m0 = _mm_loadu_ps(&in[8 * i + 0]); + m1 = _mm_loadu_ps(&in[8 * i + 4]); + m2 = _mm_load1_ps(&scale); + + /* Scale */ + m0 = _mm_mul_ps(m0, m2); + m1 = _mm_mul_ps(m1, m2); + + /* Convert */ + m4 = _mm_cvtps_epi32(m0); + m5 = _mm_cvtps_epi32(m1); + + /* Pack and store */ + m5 = _mm_packs_epi32(m4, m5); + _mm_storeu_si128((__m128i *) & out[8 * i], m5); + } +} + +/* 8*N single precision floats scaled and converted with remainder */ +void _sse_convert_scale_ps_si16(short *restrict out, + const float *restrict in, float scale, int len) +{ + int start = len / 8 * 8; + + _sse_convert_scale_ps_si16_8n(out, in, scale, len); + + for (int i = 0; i < len % 8; i++) + out[start + i] = in[start + i] * scale; +} + +/* 16*N single precision floats scaled and converted to 16-bit signed integer */ +void _sse_convert_scale_ps_si16_16n(short *restrict out, + const float *restrict in, + float scale, int len) +{ + __m128 m0, m1, m2, m3, m4; + __m128i m5, m6, m7, m8; + + for (int i = 0; i < len / 16; i++) { + /* Load (unaligned) packed floats */ + m0 = _mm_loadu_ps(&in[16 * i + 0]); + m1 = _mm_loadu_ps(&in[16 * i + 4]); + m2 = _mm_loadu_ps(&in[16 * i + 8]); + m3 = _mm_loadu_ps(&in[16 * i + 12]); + m4 = _mm_load1_ps(&scale); + + /* Scale */ + m0 = _mm_mul_ps(m0, m4); + m1 = _mm_mul_ps(m1, m4); + m2 = _mm_mul_ps(m2, m4); + m3 = _mm_mul_ps(m3, m4); + + /* Convert */ + m5 = _mm_cvtps_epi32(m0); + m6 = _mm_cvtps_epi32(m1); + m7 = _mm_cvtps_epi32(m2); + m8 = _mm_cvtps_epi32(m3); + + /* Pack and store */ + m5 = _mm_packs_epi32(m5, m6); + m7 = _mm_packs_epi32(m7, m8); + _mm_storeu_si128((__m128i *) & out[16 * i + 0], m5); + _mm_storeu_si128((__m128i *) & out[16 * i + 8], m7); + } +} +#endif diff --git a/Transceiver52M/arch/x86/convert_sse_3.h b/Transceiver52M/arch/x86/convert_sse_3.h new file mode 100644 index 0000000..c2f87d7 --- /dev/null +++ b/Transceiver52M/arch/x86/convert_sse_3.h @@ -0,0 +1,34 @@ +/* + * SSE type conversions + * Copyright (C) 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#pragma once + +/* 8*N single precision floats scaled and converted to 16-bit signed integer */ +void _sse_convert_scale_ps_si16_8n(short *restrict out, + const float *restrict in, + float scale, int len); + +/* 8*N single precision floats scaled and converted with remainder */ +void _sse_convert_scale_ps_si16(short *restrict out, + const float *restrict in, float scale, int len); + +/* 16*N single precision floats scaled and converted to 16-bit signed integer */ +void _sse_convert_scale_ps_si16_16n(short *restrict out, + const float *restrict in, + float scale, int len); diff --git a/Transceiver52M/arch/x86/convert_sse_4_1.c b/Transceiver52M/arch/x86/convert_sse_4_1.c new file mode 100644 index 0000000..42a235c --- /dev/null +++ b/Transceiver52M/arch/x86/convert_sse_4_1.c @@ -0,0 +1,77 @@ +/* + * SSE type conversions + * Copyright (C) 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include "convert_sse_4_1.h" + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_SSE4_1 +#include <smmintrin.h> + +/* 16*N 16-bit signed integer converted to single precision floats */ +void _sse_convert_si16_ps_16n(float *restrict out, + const short *restrict in, int len) +{ + __m128i m0, m1, m2, m3, m4, m5; + __m128 m6, m7, m8, m9; + + for (int i = 0; i < len / 16; i++) { + /* Load (unaligned) packed floats */ + m0 = _mm_loadu_si128((__m128i *) & in[16 * i + 0]); + m1 = _mm_loadu_si128((__m128i *) & in[16 * i + 8]); + + /* Unpack */ + m2 = _mm_cvtepi16_epi32(m0); + m4 = _mm_cvtepi16_epi32(m1); + m0 = _mm_shuffle_epi32(m0, _MM_SHUFFLE(1, 0, 3, 2)); + m1 = _mm_shuffle_epi32(m1, _MM_SHUFFLE(1, 0, 3, 2)); + m3 = _mm_cvtepi16_epi32(m0); + m5 = _mm_cvtepi16_epi32(m1); + + /* Convert */ + m6 = _mm_cvtepi32_ps(m2); + m7 = _mm_cvtepi32_ps(m3); + m8 = _mm_cvtepi32_ps(m4); + m9 = _mm_cvtepi32_ps(m5); + + /* Store */ + _mm_storeu_ps(&out[16 * i + 0], m6); + _mm_storeu_ps(&out[16 * i + 4], m7); + _mm_storeu_ps(&out[16 * i + 8], m8); + _mm_storeu_ps(&out[16 * i + 12], m9); + } +} + +/* 16*N 16-bit signed integer conversion with remainder */ +void _sse_convert_si16_ps(float *restrict out, + const short *restrict in, int len) +{ + int start = len / 16 * 16; + + _sse_convert_si16_ps_16n(out, in, len); + + for (int i = 0; i < len % 16; i++) + out[start + i] = in[start + i]; +} + +#endif diff --git a/Transceiver52M/arch/x86/convert_sse_4_1.h b/Transceiver52M/arch/x86/convert_sse_4_1.h new file mode 100644 index 0000000..57a5efb --- /dev/null +++ b/Transceiver52M/arch/x86/convert_sse_4_1.h @@ -0,0 +1,28 @@ +/* + * SSE type conversions + * Copyright (C) 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#pragma once + +/* 16*N 16-bit signed integer converted to single precision floats */ +void _sse_convert_si16_ps_16n(float *restrict out, + const short *restrict in, int len); + +/* 16*N 16-bit signed integer conversion with remainder */ +void _sse_convert_si16_ps(float *restrict out, + const short *restrict in, int len); diff --git a/Transceiver52M/arch/x86/convolve.c b/Transceiver52M/arch/x86/convolve.c new file mode 100644 index 0000000..eb38f64 --- /dev/null +++ b/Transceiver52M/arch/x86/convolve.c @@ -0,0 +1,172 @@ +/* + * SSE Convolution + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include <stdio.h> +#include "convolve.h" +#include "convolve_sse_3.h" + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +/* Architecture dependant function pointers */ +struct convolve_cpu_context { + void (*conv_cmplx_4n) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_cmplx_8n) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_cmplx) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_real4) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_real8) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_real12) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_real16) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_real20) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_real4n) (const float *, int, const float *, int, float *, + int, int, int, int, int); + void (*conv_real) (const float *, int, const float *, int, float *, int, + int, int, int, int); +}; +static struct convolve_cpu_context c; + +/* Forward declarations from base implementation */ +int _base_convolve_real(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +int _base_convolve_complex(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, + int step, int offset); + +int bounds_check(int x_len, int h_len, int y_len, + int start, int len, int step); + +/* API: Initalize convolve module */ +void convolve_init(void) +{ + c.conv_cmplx_4n = (void *)_base_convolve_complex; + c.conv_cmplx_8n = (void *)_base_convolve_complex; + c.conv_cmplx = (void *)_base_convolve_complex; + c.conv_real4 = (void *)_base_convolve_real; + c.conv_real8 = (void *)_base_convolve_real; + c.conv_real12 = (void *)_base_convolve_real; + c.conv_real16 = (void *)_base_convolve_real; + c.conv_real20 = (void *)_base_convolve_real; + c.conv_real4n = (void *)_base_convolve_real; + c.conv_real = (void *)_base_convolve_real; + +#if defined(HAVE_SSE3) && defined(HAVE___BUILTIN_CPU_SUPPORTS) + if (__builtin_cpu_supports("sse3")) { + c.conv_cmplx_4n = sse_conv_cmplx_4n; + c.conv_cmplx_8n = sse_conv_cmplx_8n; + c.conv_real4 = sse_conv_real4; + c.conv_real8 = sse_conv_real8; + c.conv_real12 = sse_conv_real12; + c.conv_real16 = sse_conv_real16; + c.conv_real20 = sse_conv_real20; + c.conv_real4n = sse_conv_real4n; + } +#endif +} + +/* API: Aligned complex-real */ +int convolve_real(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, int start, int len, int step, int offset) +{ + if (bounds_check(x_len, h_len, y_len, start, len, step) < 0) + return -1; + + memset(y, 0, len * 2 * sizeof(float)); + + if (step <= 4) { + switch (h_len) { + case 4: + c.conv_real4(x, x_len, h, h_len, y, y_len, start, len, + step, offset); + break; + case 8: + c.conv_real8(x, x_len, h, h_len, y, y_len, start, len, + step, offset); + break; + case 12: + c.conv_real12(x, x_len, h, h_len, y, y_len, start, len, + step, offset); + break; + case 16: + c.conv_real16(x, x_len, h, h_len, y, y_len, start, len, + step, offset); + break; + case 20: + c.conv_real20(x, x_len, h, h_len, y, y_len, start, len, + step, offset); + break; + default: + if (!(h_len % 4)) + c.conv_real4n(x, x_len, h, h_len, y, y_len, + start, len, step, offset); + else + c.conv_real(x, x_len, h, h_len, y, y_len, start, + len, step, offset); + } + } else + c.conv_real(x, x_len, h, h_len, y, y_len, start, len, step, + offset); + + return len; +} + +/* API: Aligned complex-complex */ +int convolve_complex(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + if (bounds_check(x_len, h_len, y_len, start, len, step) < 0) + return -1; + + memset(y, 0, len * 2 * sizeof(float)); + + if (step <= 4) { + if (!(h_len % 8)) + c.conv_cmplx_8n(x, x_len, h, h_len, y, y_len, start, + len, step, offset); + else if (!(h_len % 4)) + c.conv_cmplx_4n(x, x_len, h, h_len, y, y_len, start, + len, step, offset); + else + c.conv_cmplx(x, x_len, h, h_len, y, y_len, start, len, + step, offset); + } else + c.conv_cmplx(x, x_len, h, h_len, y, y_len, start, len, step, + offset); + + return len; +} diff --git a/Transceiver52M/arch/x86/convolve_sse_3.c b/Transceiver52M/arch/x86/convolve_sse_3.c new file mode 100644 index 0000000..dbee3cc --- /dev/null +++ b/Transceiver52M/arch/x86/convolve_sse_3.c @@ -0,0 +1,542 @@ +/* + * SSE Convolution + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#include <malloc.h> +#include <string.h> +#include <stdio.h> +#include "convolve_sse_3.h" + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#ifdef HAVE_SSE3 +#include <xmmintrin.h> +#include <pmmintrin.h> + +/* 4-tap SSE complex-real convolution */ +void sse_conv_real4(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* NOTE: The parameter list of this function has to match the parameter + * list of _base_convolve_real() in convolve_base.c. This specific + * implementation, ignores some of the parameters of + * _base_convolve_complex(), which are: x_len, y_len, offset, step */ + + __m128 m0, m1, m2, m3, m4, m5, m6, m7; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[0]); + m1 = _mm_load_ps(&h[4]); + m7 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + + for (int i = 0; i < len; i++) { + /* Load (unaligned) input data */ + m0 = _mm_loadu_ps(&_x[2 * i + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 4]); + m2 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m3 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Quad multiply */ + m4 = _mm_mul_ps(m2, m7); + m5 = _mm_mul_ps(m3, m7); + + /* Sum and store */ + m6 = _mm_hadd_ps(m4, m5); + m0 = _mm_hadd_ps(m6, m6); + + _mm_store_ss(&y[2 * i + 0], m0); + m0 = _mm_shuffle_ps(m0, m0, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m0); + } +} + +/* 8-tap SSE complex-real convolution */ +void sse_conv_real8(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* See NOTE in sse_conv_real4() */ + + __m128 m0, m1, m2, m3, m4, m5, m6, m7, m8, m9; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[0]); + m1 = _mm_load_ps(&h[4]); + m2 = _mm_load_ps(&h[8]); + m3 = _mm_load_ps(&h[12]); + + m4 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m5 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + + for (int i = 0; i < len; i++) { + /* Load (unaligned) input data */ + m0 = _mm_loadu_ps(&_x[2 * i + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 4]); + m2 = _mm_loadu_ps(&_x[2 * i + 8]); + m3 = _mm_loadu_ps(&_x[2 * i + 12]); + + m6 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m7 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m8 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m9 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Quad multiply */ + m6 = _mm_mul_ps(m6, m4); + m7 = _mm_mul_ps(m7, m4); + m8 = _mm_mul_ps(m8, m5); + m9 = _mm_mul_ps(m9, m5); + + /* Sum and store */ + m6 = _mm_add_ps(m6, m8); + m7 = _mm_add_ps(m7, m9); + m6 = _mm_hadd_ps(m6, m7); + m6 = _mm_hadd_ps(m6, m6); + + _mm_store_ss(&y[2 * i + 0], m6); + m6 = _mm_shuffle_ps(m6, m6, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m6); + } +} + +/* 12-tap SSE complex-real convolution */ +void sse_conv_real12(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* See NOTE in sse_conv_real4() */ + + __m128 m0, m1, m2, m3, m4, m5, m6, m7; + __m128 m8, m9, m10, m11, m12, m13, m14; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[0]); + m1 = _mm_load_ps(&h[4]); + m2 = _mm_load_ps(&h[8]); + m3 = _mm_load_ps(&h[12]); + m4 = _mm_load_ps(&h[16]); + m5 = _mm_load_ps(&h[20]); + + m12 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m13 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m14 = _mm_shuffle_ps(m4, m5, _MM_SHUFFLE(0, 2, 0, 2)); + + for (int i = 0; i < len; i++) { + /* Load (unaligned) input data */ + m0 = _mm_loadu_ps(&_x[2 * i + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 4]); + m2 = _mm_loadu_ps(&_x[2 * i + 8]); + m3 = _mm_loadu_ps(&_x[2 * i + 12]); + + m4 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m5 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m6 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m7 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + + m0 = _mm_loadu_ps(&_x[2 * i + 16]); + m1 = _mm_loadu_ps(&_x[2 * i + 20]); + + m8 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m9 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Quad multiply */ + m0 = _mm_mul_ps(m4, m12); + m1 = _mm_mul_ps(m5, m12); + m2 = _mm_mul_ps(m6, m13); + m3 = _mm_mul_ps(m7, m13); + m4 = _mm_mul_ps(m8, m14); + m5 = _mm_mul_ps(m9, m14); + + /* Sum and store */ + m8 = _mm_add_ps(m0, m2); + m9 = _mm_add_ps(m1, m3); + m10 = _mm_add_ps(m8, m4); + m11 = _mm_add_ps(m9, m5); + + m2 = _mm_hadd_ps(m10, m11); + m3 = _mm_hadd_ps(m2, m2); + + _mm_store_ss(&y[2 * i + 0], m3); + m3 = _mm_shuffle_ps(m3, m3, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m3); + } +} + +/* 16-tap SSE complex-real convolution */ +void sse_conv_real16(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* See NOTE in sse_conv_real4() */ + + __m128 m0, m1, m2, m3, m4, m5, m6, m7; + __m128 m8, m9, m10, m11, m12, m13, m14, m15; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[0]); + m1 = _mm_load_ps(&h[4]); + m2 = _mm_load_ps(&h[8]); + m3 = _mm_load_ps(&h[12]); + + m4 = _mm_load_ps(&h[16]); + m5 = _mm_load_ps(&h[20]); + m6 = _mm_load_ps(&h[24]); + m7 = _mm_load_ps(&h[28]); + + m12 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m13 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m14 = _mm_shuffle_ps(m4, m5, _MM_SHUFFLE(0, 2, 0, 2)); + m15 = _mm_shuffle_ps(m6, m7, _MM_SHUFFLE(0, 2, 0, 2)); + + for (int i = 0; i < len; i++) { + /* Load (unaligned) input data */ + m0 = _mm_loadu_ps(&_x[2 * i + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 4]); + m2 = _mm_loadu_ps(&_x[2 * i + 8]); + m3 = _mm_loadu_ps(&_x[2 * i + 12]); + + m4 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m5 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m6 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m7 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + + m0 = _mm_loadu_ps(&_x[2 * i + 16]); + m1 = _mm_loadu_ps(&_x[2 * i + 20]); + m2 = _mm_loadu_ps(&_x[2 * i + 24]); + m3 = _mm_loadu_ps(&_x[2 * i + 28]); + + m8 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m9 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m10 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m11 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Quad multiply */ + m0 = _mm_mul_ps(m4, m12); + m1 = _mm_mul_ps(m5, m12); + m2 = _mm_mul_ps(m6, m13); + m3 = _mm_mul_ps(m7, m13); + + m4 = _mm_mul_ps(m8, m14); + m5 = _mm_mul_ps(m9, m14); + m6 = _mm_mul_ps(m10, m15); + m7 = _mm_mul_ps(m11, m15); + + /* Sum and store */ + m8 = _mm_add_ps(m0, m2); + m9 = _mm_add_ps(m1, m3); + m10 = _mm_add_ps(m4, m6); + m11 = _mm_add_ps(m5, m7); + + m0 = _mm_add_ps(m8, m10); + m1 = _mm_add_ps(m9, m11); + m2 = _mm_hadd_ps(m0, m1); + m3 = _mm_hadd_ps(m2, m2); + + _mm_store_ss(&y[2 * i + 0], m3); + m3 = _mm_shuffle_ps(m3, m3, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m3); + } +} + +/* 20-tap SSE complex-real convolution */ +void sse_conv_real20(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* See NOTE in sse_conv_real4() */ + + __m128 m0, m1, m2, m3, m4, m5, m6, m7; + __m128 m8, m9, m11, m12, m13, m14, m15; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[0]); + m1 = _mm_load_ps(&h[4]); + m2 = _mm_load_ps(&h[8]); + m3 = _mm_load_ps(&h[12]); + m4 = _mm_load_ps(&h[16]); + m5 = _mm_load_ps(&h[20]); + m6 = _mm_load_ps(&h[24]); + m7 = _mm_load_ps(&h[28]); + m8 = _mm_load_ps(&h[32]); + m9 = _mm_load_ps(&h[36]); + + m11 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m12 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m13 = _mm_shuffle_ps(m4, m5, _MM_SHUFFLE(0, 2, 0, 2)); + m14 = _mm_shuffle_ps(m6, m7, _MM_SHUFFLE(0, 2, 0, 2)); + m15 = _mm_shuffle_ps(m8, m9, _MM_SHUFFLE(0, 2, 0, 2)); + + for (int i = 0; i < len; i++) { + /* Multiply-accumulate first 12 taps */ + m0 = _mm_loadu_ps(&_x[2 * i + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 4]); + m2 = _mm_loadu_ps(&_x[2 * i + 8]); + m3 = _mm_loadu_ps(&_x[2 * i + 12]); + m4 = _mm_loadu_ps(&_x[2 * i + 16]); + m5 = _mm_loadu_ps(&_x[2 * i + 20]); + + m6 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m7 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m8 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m9 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + m0 = _mm_shuffle_ps(m4, m5, _MM_SHUFFLE(0, 2, 0, 2)); + m1 = _mm_shuffle_ps(m4, m5, _MM_SHUFFLE(1, 3, 1, 3)); + + m2 = _mm_mul_ps(m6, m11); + m3 = _mm_mul_ps(m7, m11); + m4 = _mm_mul_ps(m8, m12); + m5 = _mm_mul_ps(m9, m12); + m6 = _mm_mul_ps(m0, m13); + m7 = _mm_mul_ps(m1, m13); + + m0 = _mm_add_ps(m2, m4); + m1 = _mm_add_ps(m3, m5); + m8 = _mm_add_ps(m0, m6); + m9 = _mm_add_ps(m1, m7); + + /* Multiply-accumulate last 8 taps */ + m0 = _mm_loadu_ps(&_x[2 * i + 24]); + m1 = _mm_loadu_ps(&_x[2 * i + 28]); + m2 = _mm_loadu_ps(&_x[2 * i + 32]); + m3 = _mm_loadu_ps(&_x[2 * i + 36]); + + m4 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m5 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m6 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m7 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + + m0 = _mm_mul_ps(m4, m14); + m1 = _mm_mul_ps(m5, m14); + m2 = _mm_mul_ps(m6, m15); + m3 = _mm_mul_ps(m7, m15); + + m4 = _mm_add_ps(m0, m2); + m5 = _mm_add_ps(m1, m3); + + /* Final sum and store */ + m0 = _mm_add_ps(m8, m4); + m1 = _mm_add_ps(m9, m5); + m2 = _mm_hadd_ps(m0, m1); + m3 = _mm_hadd_ps(m2, m2); + + _mm_store_ss(&y[2 * i + 0], m3); + m3 = _mm_shuffle_ps(m3, m3, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m3); + } +} + +/* 4*N-tap SSE complex-real convolution */ +void sse_conv_real4n(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* See NOTE in sse_conv_real4() */ + + __m128 m0, m1, m2, m4, m5, m6, m7; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + for (int i = 0; i < len; i++) { + /* Zero */ + m6 = _mm_setzero_ps(); + m7 = _mm_setzero_ps(); + + for (int n = 0; n < h_len / 4; n++) { + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[8 * n + 0]); + m1 = _mm_load_ps(&h[8 * n + 4]); + m2 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + + /* Load (unaligned) input data */ + m0 = _mm_loadu_ps(&_x[2 * i + 8 * n + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 8 * n + 4]); + m4 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m5 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Quad multiply */ + m0 = _mm_mul_ps(m2, m4); + m1 = _mm_mul_ps(m2, m5); + + /* Accumulate */ + m6 = _mm_add_ps(m6, m0); + m7 = _mm_add_ps(m7, m1); + } + + m0 = _mm_hadd_ps(m6, m7); + m0 = _mm_hadd_ps(m0, m0); + + _mm_store_ss(&y[2 * i + 0], m0); + m0 = _mm_shuffle_ps(m0, m0, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m0); + } +} + +/* 4*N-tap SSE complex-complex convolution */ +void sse_conv_cmplx_4n(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* NOTE: The parameter list of this function has to match the parameter + * list of _base_convolve_complex() in convolve_base.c. This specific + * implementation, ignores some of the parameters of + * _base_convolve_complex(), which are: x_len, y_len, offset, step. */ + + __m128 m0, m1, m2, m3, m4, m5, m6, m7; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + for (int i = 0; i < len; i++) { + /* Zero */ + m6 = _mm_setzero_ps(); + m7 = _mm_setzero_ps(); + + for (int n = 0; n < h_len / 4; n++) { + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[8 * n + 0]); + m1 = _mm_load_ps(&h[8 * n + 4]); + m2 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m3 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Load (unaligned) input data */ + m0 = _mm_loadu_ps(&_x[2 * i + 8 * n + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 8 * n + 4]); + m4 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m5 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Quad multiply */ + m0 = _mm_mul_ps(m2, m4); + m1 = _mm_mul_ps(m3, m5); + + m2 = _mm_mul_ps(m2, m5); + m3 = _mm_mul_ps(m3, m4); + + /* Sum */ + m0 = _mm_sub_ps(m0, m1); + m2 = _mm_add_ps(m2, m3); + + /* Accumulate */ + m6 = _mm_add_ps(m6, m0); + m7 = _mm_add_ps(m7, m2); + } + + m0 = _mm_hadd_ps(m6, m7); + m0 = _mm_hadd_ps(m0, m0); + + _mm_store_ss(&y[2 * i + 0], m0); + m0 = _mm_shuffle_ps(m0, m0, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m0); + } +} + +/* 8*N-tap SSE complex-complex convolution */ +void sse_conv_cmplx_8n(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset) +{ + /* See NOTE in sse_conv_cmplx_4n() */ + + __m128 m0, m1, m2, m3, m4, m5, m6, m7; + __m128 m8, m9, m10, m11, m12, m13, m14, m15; + + const float *_x = &x[2 * (-(h_len - 1) + start)]; + + for (int i = 0; i < len; i++) { + /* Zero */ + m12 = _mm_setzero_ps(); + m13 = _mm_setzero_ps(); + m14 = _mm_setzero_ps(); + m15 = _mm_setzero_ps(); + + for (int n = 0; n < h_len / 8; n++) { + /* Load (aligned) filter taps */ + m0 = _mm_load_ps(&h[16 * n + 0]); + m1 = _mm_load_ps(&h[16 * n + 4]); + m2 = _mm_load_ps(&h[16 * n + 8]); + m3 = _mm_load_ps(&h[16 * n + 12]); + + m4 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m5 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m6 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m7 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Load (unaligned) input data */ + m0 = _mm_loadu_ps(&_x[2 * i + 16 * n + 0]); + m1 = _mm_loadu_ps(&_x[2 * i + 16 * n + 4]); + m2 = _mm_loadu_ps(&_x[2 * i + 16 * n + 8]); + m3 = _mm_loadu_ps(&_x[2 * i + 16 * n + 12]); + + m8 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(0, 2, 0, 2)); + m9 = _mm_shuffle_ps(m0, m1, _MM_SHUFFLE(1, 3, 1, 3)); + m10 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(0, 2, 0, 2)); + m11 = _mm_shuffle_ps(m2, m3, _MM_SHUFFLE(1, 3, 1, 3)); + + /* Quad multiply */ + m0 = _mm_mul_ps(m4, m8); + m1 = _mm_mul_ps(m5, m9); + m2 = _mm_mul_ps(m6, m10); + m3 = _mm_mul_ps(m7, m11); + + m4 = _mm_mul_ps(m4, m9); + m5 = _mm_mul_ps(m5, m8); + m6 = _mm_mul_ps(m6, m11); + m7 = _mm_mul_ps(m7, m10); + + /* Sum */ + m0 = _mm_sub_ps(m0, m1); + m2 = _mm_sub_ps(m2, m3); + m4 = _mm_add_ps(m4, m5); + m6 = _mm_add_ps(m6, m7); + + /* Accumulate */ + m12 = _mm_add_ps(m12, m0); + m13 = _mm_add_ps(m13, m2); + m14 = _mm_add_ps(m14, m4); + m15 = _mm_add_ps(m15, m6); + } + + m0 = _mm_add_ps(m12, m13); + m1 = _mm_add_ps(m14, m15); + m2 = _mm_hadd_ps(m0, m1); + m2 = _mm_hadd_ps(m2, m2); + + _mm_store_ss(&y[2 * i + 0], m2); + m2 = _mm_shuffle_ps(m2, m2, _MM_SHUFFLE(0, 3, 2, 1)); + _mm_store_ss(&y[2 * i + 1], m2); + } +} +#endif diff --git a/Transceiver52M/arch/x86/convolve_sse_3.h b/Transceiver52M/arch/x86/convolve_sse_3.h new file mode 100644 index 0000000..ac30ca5 --- /dev/null +++ b/Transceiver52M/arch/x86/convolve_sse_3.h @@ -0,0 +1,68 @@ +/* + * SSE Convolution + * Copyright (C) 2012, 2013 Thomas Tsou <tom@tsou.cc> + * + * 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 + */ + +#pragma once + +/* 4-tap SSE complex-real convolution */ +void sse_conv_real4(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); + +/* 8-tap SSE complex-real convolution */ +void sse_conv_real8(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); + +/* 12-tap SSE complex-real convolution */ +void sse_conv_real12(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); + +/* 16-tap SSE complex-real convolution */ +void sse_conv_real16(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); + +/* 20-tap SSE complex-real convolution */ +void sse_conv_real20(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); + +/* 4*N-tap SSE complex-real convolution */ +void sse_conv_real4n(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); + +/* 4*N-tap SSE complex-complex convolution */ +void sse_conv_cmplx_4n(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); + +/* 8*N-tap SSE complex-complex convolution */ +void sse_conv_cmplx_8n(const float *x, int x_len, + const float *h, int h_len, + float *y, int y_len, + int start, int len, int step, int offset); |