git.sesse.net Git - vlc/blob - modules/codec/a52old/imdct/srfft_3dn.c

   1 /*****************************************************************************
   2  * srfft_3dn.c: accelerated 3D Now! A52 fft functions
   3  *****************************************************************************
   4  * Copyright (C) 1999-2001 VideoLAN
   5  * $Id: srfft_3dn.c,v 1.1 2002/08/04 17:23:42 sam Exp $
   6  *
   7  * Authors: Renaud Dartus <reno@videolan.org>
   8  *
   9  * This program is free software; you can redistribute it and/or modify
  10  * it under the terms of the GNU General Public License as published by
  11  * the Free Software Foundation; either version 2 of the License, or
  12  * (at your option) any later version.
  13  *
  14  * This program is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  * GNU General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU General Public License
  20  * along with this program; if not, write to the Free Software
  21  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111, USA.
  22  *****************************************************************************/
  23
  24 /*****************************************************************************
  25  * Preamble
  26  *****************************************************************************/
  27 #include <math.h>
  28 #include <stdio.h>
  29
  30 #include <vlc/vlc.h>
  31
  32 #include "../imdct.h"
  33 #include "srfft.h"
  34
  35 void hsqrt2_3dn (void);
  36 void C_1_3dn (void);
  37 static void fft_4_3dn (complex_t *x);
  38 static void fft_8_3dn (complex_t *x);
  39 static void fft_asmb_3dn (int k, complex_t *x, complex_t *wTB,
  40              const complex_t *d, const complex_t *d_3);
  41
  42 void E_( fft_64p ) ( complex_t *a )
  43 {
  44         fft_8_3dn(&a[0]); fft_4_3dn(&a[8]); fft_4_3dn(&a[12]);
  45         fft_asmb_3dn(2, &a[0], &a[8], &delta16[0], &delta16_3[0]);
  46
  47         fft_8_3dn(&a[16]), fft_8_3dn(&a[24]);
  48         fft_asmb_3dn(4, &a[0], &a[16],&delta32[0], &delta32_3[0]);
  49
  50         fft_8_3dn(&a[32]); fft_4_3dn(&a[40]); fft_4_3dn(&a[44]);
  51         fft_asmb_3dn(2, &a[32], &a[40], &delta16[0], &delta16_3[0]);
  52
  53         fft_8_3dn(&a[48]); fft_4_3dn(&a[56]); fft_4_3dn(&a[60]);
  54         fft_asmb_3dn(2, &a[48], &a[56], &delta16[0], &delta16_3[0]);
  55
  56         fft_asmb_3dn(8, &a[0], &a[32],&delta64[0], &delta64_3[0]);
  57 }
  58
  59 void E_( fft_128p ) ( complex_t *a )
  60 {
  61     fft_8_3dn(&a[0]); fft_4_3dn(&a[8]); fft_4_3dn(&a[12]);
  62         fft_asmb_3dn(2, &a[0], &a[8], &delta16[0], &delta16_3[0]);
  63
  64         fft_8_3dn(&a[16]), fft_8_3dn(&a[24]);
  65         fft_asmb_3dn(4, &a[0], &a[16],&delta32[0], &delta32_3[0]);
  66
  67         fft_8_3dn(&a[32]); fft_4_3dn(&a[40]); fft_4_3dn(&a[44]);
  68         fft_asmb_3dn(2, &a[32], &a[40], &delta16[0], &delta16_3[0]);
  69
  70         fft_8_3dn(&a[48]); fft_4_3dn(&a[56]); fft_4_3dn(&a[60]);
  71         fft_asmb_3dn(2, &a[48], &a[56], &delta16[0], &delta16_3[0]);
  72
  73         fft_asmb_3dn(8, &a[0], &a[32],&delta64[0], &delta64_3[0]);
  74
  75         fft_8_3dn(&a[64]); fft_4_3dn(&a[72]); fft_4_3dn(&a[76]);
  76         /* fft_16(&a[64]); */
  77         fft_asmb_3dn(2, &a[64], &a[72], &delta16[0], &delta16_3[0]);
  78
  79         fft_8_3dn(&a[80]); fft_8_3dn(&a[88]);
  80
  81         /* fft_32(&a[64]); */
  82         fft_asmb_3dn(4, &a[64], &a[80],&delta32[0], &delta32_3[0]);
  83
  84         fft_8_3dn(&a[96]); fft_4_3dn(&a[104]), fft_4_3dn(&a[108]);
  85         /* fft_16(&a[96]); */
  86         fft_asmb_3dn(2, &a[96], &a[104], &delta16[0], &delta16_3[0]);
  87
  88         fft_8_3dn(&a[112]), fft_8_3dn(&a[120]);
  89         /* fft_32(&a[96]); */
  90         fft_asmb_3dn(4, &a[96], &a[112], &delta32[0], &delta32_3[0]);
  91
  92         /* fft_128(&a[0]); */
  93         fft_asmb_3dn(16, &a[0], &a[64], &delta128[0], &delta128_3[0]);
  94 }
  95
  96 void hsqrt2_3dn (void)
  97 {
  98     __asm__ (
  99      ".float 0f0.707106781188\n"
 100      ".float 0f0.707106781188\n"
 101      ".float 0f-0.707106781188\n"
 102      ".float 0f-0.707106781188\n"
 103      );
 104 }
 105
 106 void C_1_3dn (void)
 107 {
 108     __asm__ (
 109      ".float 0f-1.0\n"
 110      ".float 0f1.0\n"
 111      ".float 0f-1.0\n"
 112      ".float 0f1.0\n"
 113      );
 114 }
 115
 116 static void fft_4_3dn (complex_t *x)
 117 {
 118     __asm__ __volatile__ (
 119     ".align 16\n"
 120         "movq    (%%eax), %%mm0\n"      /* x[0] */
 121         "movq   8(%%eax), %%mm1\n"      /* x[1] */
 122         "movq  16(%%eax), %%mm2\n"      /* x[2] */
 123         "movq  24(%%eax), %%mm3\n"      /* x[3] */
 124         "movq    %%mm0, %%mm4\n"            /* x[1] */
 125         "movq    %%mm1, %%mm5\n"                /* x[1] */
 126         "movq    %%mm0, %%mm6\n"                /* x[0] */
 127         "pfadd   %%mm2, %%mm0\n"                /* x[0] + x[2] */
 128         "pfadd   %%mm3, %%mm1\n"                /* x[1] + x[3] */
 129         "pfsub   %%mm2, %%mm4\n"                /* x[0] - x[2] */
 130         "pfsub   %%mm3, %%mm5\n"                /* x[1] - x[3] */
 131
 132     "pfadd   %%mm1, %%mm0\n"        /* x[0] + x[2] + x[1] + x[3] */
 133     "pfsub   %%mm1, %%mm6\n"        /* x[0] + x[2] - x[1] - x[3] */
 134
 135     "movq   %%mm0, (%%eax)\n"
 136     "movq   %%mm6, 16(%%eax)\n"
 137
 138     "pxor    %%mm6, %%mm6\n"
 139     "movq    %%mm5, %%mm2\n"        /* x[1] - x[3] */
 140     "movq    %%mm4, %%mm3\n"        /* x[0] - x[2] */
 141     "pfsub   %%mm5, %%mm6\n"        /* x[3] - x[1] */
 142
 143     "punpckhdq %%mm2,%%mm2\n"       /* x[1] - x[3].im */
 144     "punpckldq %%mm6,%%mm6\n"       /* x[3] - x[1].re */
 145     "punpckhdq %%mm6,%%mm2\n"       /* x[3] - x[1].re,  x[1] - x[3].im */
 146
 147         "pfsub   %%mm2, %%mm4\n"        /* x0i-x2i-x3r+x1.r,x0r-x2r-x1i+x3i */
 148     "pfadd   %%mm3, %%mm2\n"        /* x0i-x2i+x3r-x1.r, x0r-x2r+x1i-x3.i */
 149
 150     "movq  %%mm2,  8(%%eax)\n"    /* mm4_2 + mm6_1, mm4_1 + mm5_2 */
 151         "movq  %%mm4, 24(%%eax)\n"    /* mm4_2 - mm6_1, mm4_1 - mm5_2 */
 152         "femms\n"
 153     : "=a" (x)
 154     : "a" (x) );
 155 }
 156
 157 static void fft_8_3dn (complex_t *x)
 158 {
 159   register float wT1_r, wT1_i, wB1_r, wB1_i, wT2_r, wT2_i, wB2_r, wB2_i;
 160
 161   wT1_r = x[1].real;
 162   wT1_i = x[1].imag;
 163   wB1_r = x[3].real;
 164   wB1_i = x[3].imag;
 165
 166   x[1] = x[2];
 167   x[2] = x[4];
 168   x[3] = x[6];
 169   { /* fft_4 */
 170       register float yt_r, yt_i, yb_r, yb_i, u_r, u_i, vi_r, vi_i;
 171
 172       yt_r = x[0].real;
 173       yb_r = yt_r - x[2].real;
 174       yt_r += x[2].real;
 175
 176       u_r = x[1].real;
 177       vi_i = x[3].real - u_r;
 178       u_r += x[3].real;
 179
 180       u_i = x[1].imag;
 181       vi_r = u_i - x[3].imag;
 182       u_i += x[3].imag;
 183
 184       yt_i = yt_r;
 185       yt_i += u_r;
 186       x[0].real = yt_i;
 187       yt_r -= u_r;
 188       x[2].real = yt_r;
 189       yt_i = yb_r;
 190       yt_i += vi_r;
 191       x[1].real = yt_i;
 192       yb_r -= vi_r;
 193       x[3].real = yb_r;
 194
 195       yt_i = x[0].imag;
 196       yb_i = yt_i - x[2].imag;
 197       yt_i += x[2].imag;
 198
 199       yt_r = yt_i;
 200       yt_r += u_i;
 201       x[0].imag = yt_r;
 202       yt_i -= u_i;
 203       x[2].imag = yt_i;
 204       yt_r = yb_i;
 205       yt_r += vi_i;
 206       x[1].imag = yt_r;
 207       yb_i -= vi_i;
 208       x[3].imag = yb_i;
 209   }
 210
 211   /* x[0] x[4] */
 212   wT2_r = x[5].real;
 213   wT2_r += x[7].real;
 214   wT2_r += wT1_r;
 215   wT2_r += wB1_r;
 216   wT2_i = wT2_r;
 217   wT2_r += x[0].real;
 218   wT2_i = x[0].real - wT2_i;
 219   x[0].real = wT2_r;
 220   x[4].real = wT2_i;
 221
 222   wT2_i = x[5].imag;
 223   wT2_i += x[7].imag;
 224   wT2_i += wT1_i;
 225   wT2_i += wB1_i;
 226   wT2_r = wT2_i;
 227   wT2_r += x[0].imag;
 228   wT2_i = x[0].imag - wT2_i;
 229   x[0].imag = wT2_r;
 230   x[4].imag = wT2_i;
 231
 232   /* x[2] x[6] */
 233   wT2_r = x[5].imag;
 234   wT2_r -= x[7].imag;
 235   wT2_r += wT1_i;
 236   wT2_r -= wB1_i;
 237   wT2_i = wT2_r;
 238   wT2_r += x[2].real;
 239   wT2_i = x[2].real - wT2_i;
 240   x[2].real = wT2_r;
 241   x[6].real = wT2_i;
 242
 243   wT2_i = x[5].real;
 244   wT2_i -= x[7].real;
 245   wT2_i += wT1_r;
 246   wT2_i -= wB1_r;
 247   wT2_r = wT2_i;
 248   wT2_r += x[2].imag;
 249   wT2_i = x[2].imag - wT2_i;
 250   x[2].imag = wT2_i;
 251   x[6].imag = wT2_r;
 252
 253   /* x[1] x[5] */
 254   wT2_r = wT1_r;
 255   wT2_r += wB1_i;
 256   wT2_r -= x[5].real;
 257   wT2_r -= x[7].imag;
 258   wT2_i = wT1_i;
 259   wT2_i -= wB1_r;
 260   wT2_i -= x[5].imag;
 261   wT2_i += x[7].real;
 262
 263   wB2_r = wT2_r;
 264   wB2_r += wT2_i;
 265   wT2_i -= wT2_r;
 266   wB2_r *= HSQRT2;
 267   wT2_i *= HSQRT2;
 268   wT2_r = wB2_r;
 269   wB2_r += x[1].real;
 270   wT2_r =  x[1].real - wT2_r;
 271
 272   wB2_i = x[5].real;
 273   x[1].real = wB2_r;
 274   x[5].real = wT2_r;
 275
 276   wT2_r = wT2_i;
 277   wT2_r += x[1].imag;
 278   wT2_i = x[1].imag - wT2_i;
 279   wB2_r = x[5].imag;
 280   x[1].imag = wT2_r;
 281   x[5].imag = wT2_i;
 282
 283   /* x[3] x[7] */
 284   wT1_r -= wB1_i;
 285   wT1_i += wB1_r;
 286   wB1_r = wB2_i - x[7].imag;
 287   wB1_i = wB2_r + x[7].real;
 288   wT1_r -= wB1_r;
 289   wT1_i -= wB1_i;
 290   wB1_r = wT1_r + wT1_i;
 291   wB1_r *= HSQRT2;
 292   wT1_i -= wT1_r;
 293   wT1_i *= HSQRT2;
 294   wB2_r = x[3].real;
 295   wB2_i = wB2_r + wT1_i;
 296   wB2_r -= wT1_i;
 297   x[3].real = wB2_i;
 298   x[7].real = wB2_r;
 299   wB2_i = x[3].imag;
 300   wB2_r = wB2_i + wB1_r;
 301   wB2_i -= wB1_r;
 302   x[3].imag = wB2_i;
 303   x[7].imag = wB2_r;
 304 }
 305
 306
 307 static void fft_asmb_3dn (int k, complex_t *x, complex_t *wTB,
 308              const complex_t *d, const complex_t *d_3)
 309 {
 310   register complex_t  *x2k, *x3k, *x4k, *wB;
 311   register float a_r, a_i, a1_r, a1_i, u_r, u_i, v_r, v_i;
 312
 313   x2k = x + 2 * k;
 314   x3k = x2k + 2 * k;
 315   x4k = x3k + 2 * k;
 316   wB = wTB + 2 * k;
 317
 318   TRANSZERO(x[0],x2k[0],x3k[0],x4k[0]);
 319   TRANS(x[1],x2k[1],x3k[1],x4k[1],wTB[1],wB[1],d[1],d_3[1]);
 320
 321   --k;
 322   for(;;) {
 323      TRANS(x[2],x2k[2],x3k[2],x4k[2],wTB[2],wB[2],d[2],d_3[2]);
 324      TRANS(x[3],x2k[3],x3k[3],x4k[3],wTB[3],wB[3],d[3],d_3[3]);
 325      if (!--k) break;
 326      x += 2;
 327      x2k += 2;
 328      x3k += 2;
 329      x4k += 2;
 330      d += 2;
 331      d_3 += 2;
 332      wTB += 2;
 333      wB += 2;
 334   }
 335 }