"Fossies" - the Fresh Open Source Software Archive

Member "faac-1_30/libfaac/fft.c" (16 Oct 2019, 9331 Bytes) of package /linux/misc/faac-1_30.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "fft.c" see the Fossies "Dox" file reference documentation and the last Fossies "Diffs" side-by-side code changes report: 1.28_vs_1.29.9.2.

    1 /*
    2  * FAAC - Freeware Advanced Audio Coder
    3  * $Id: fft.c,v 1.12 2005/02/02 07:49:55 sur Exp $
    4  * Copyright (C) 2002 Krzysztof Nikiel
    5  *
    6  * This library 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  * This library 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 this library; if not, write to the Free Software
   18  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
   19  *
   20  */
   21 
   22 #include <math.h>
   23 #include <stdlib.h>
   24 #include <stdio.h>
   25 
   26 #include "fft.h"
   27 #include "util.h"
   28 
   29 #define MAXLOGM 9
   30 #define MAXLOGR 8
   31 
   32 #if defined DRM && !defined DRM_1024
   33 
   34 #include "kiss_fft/kiss_fft.h"
   35 #include "kiss_fft/kiss_fftr.h"
   36 
   37 static const int logm_to_nfft[] = 
   38 {
   39 /*  0       1       2       3       */
   40     0,      0,      0,      0,
   41 /*  4       5       6       7       */
   42     0,      0,      60,     0,
   43 /*  8       9                       */
   44     240,    480
   45 };
   46 
   47 void fft_initialize( FFT_Tables *fft_tables )
   48 {
   49     memset( fft_tables->cfg, 0, sizeof( fft_tables->cfg ) );
   50 }
   51 void fft_terminate( FFT_Tables *fft_tables )
   52 {
   53     unsigned int i;
   54     for ( i = 0; i < sizeof( fft_tables->cfg ) / sizeof( fft_tables->cfg[0] ); i++ )
   55     {
   56         if ( fft_tables->cfg[i][0] )
   57         {
   58             free( fft_tables->cfg[i][0] );
   59             fft_tables->cfg[i][0] = NULL;
   60         }
   61         if ( fft_tables->cfg[i][1] )
   62         {
   63             free( fft_tables->cfg[i][1] );
   64             fft_tables->cfg[i][1] = NULL;
   65         }
   66     }
   67 }
   68 
   69 void rfft( FFT_Tables *fft_tables, double *x, int logm )
   70 {
   71 /* sur: use real-only optimized FFT */
   72 
   73     int nfft = 0;
   74 
   75     kiss_fft_scalar fin[1 << MAXLOGR];
   76     kiss_fft_cpx    fout[1 << MAXLOGR];
   77 
   78     if ( logm > MAXLOGR )
   79     {
   80         fprintf(stderr, "fft size too big\n");
   81         exit(1);
   82     }
   83 
   84     nfft = logm_to_nfft[logm];
   85 
   86     if ( fft_tables->cfg[logm][0] == NULL )
   87     {
   88         if ( nfft )
   89         {
   90             fft_tables->cfg[logm][0] = kiss_fftr_alloc( nfft, 0, NULL, NULL );
   91         }
   92         else
   93         {
   94             fprintf(stderr, "bad logm = %d\n", logm);
   95             exit( 1 );
   96         }
   97     }
   98 
   99     if ( fft_tables->cfg[logm][0] )
  100     {
  101         unsigned int i;
  102         
  103         for ( i = 0; i < nfft; i++ )
  104         {
  105             fin[i] = x[i];    
  106         }
  107         
  108         kiss_fftr( (kiss_fftr_cfg)fft_tables->cfg[logm][0], fin, fout );
  109 
  110         for ( i = 0; i < nfft / 2; i++ )
  111         {
  112             x[i]            = fout[i].r;
  113             x[i + nfft / 2] = fout[i].i;
  114         }
  115     }
  116     else
  117     {
  118         fprintf( stderr, "bad config for logm = %d\n", logm);
  119         exit( 1 );
  120     }
  121 }
  122 
  123 void fft( FFT_Tables *fft_tables, double *xr, double *xi, int logm )
  124 {
  125     int nfft = 0;
  126 
  127     kiss_fft_cpx    fin[1 << MAXLOGM];
  128     kiss_fft_cpx    fout[1 << MAXLOGM];
  129 
  130     if ( logm > MAXLOGM )
  131     {
  132         fprintf(stderr, "fft size too big\n");
  133         exit(1);
  134     }
  135 
  136     nfft = logm_to_nfft[logm];
  137 
  138     if ( fft_tables->cfg[logm][0] == NULL )
  139     {
  140         if ( nfft )
  141         {
  142             fft_tables->cfg[logm][0] = kiss_fft_alloc( nfft, 0, NULL, NULL );
  143         }
  144         else
  145         {
  146             fprintf(stderr, "bad logm = %d\n", logm);
  147             exit( 1 );
  148         }
  149     }
  150 
  151     if ( fft_tables->cfg[logm][0] )
  152     {
  153         unsigned int i;
  154         
  155         for ( i = 0; i < nfft; i++ )
  156         {
  157             fin[i].r = xr[i];    
  158             fin[i].i = xi[i];
  159         }
  160         
  161         kiss_fft( (kiss_fft_cfg)fft_tables->cfg[logm][0], fin, fout );
  162 
  163         for ( i = 0; i < nfft; i++ )
  164         {
  165             xr[i]   = fout[i].r;
  166             xi[i]   = fout[i].i;
  167         }
  168     }
  169     else
  170     {
  171         fprintf( stderr, "bad config for logm = %d\n", logm);
  172         exit( 1 );
  173     }
  174 }
  175 
  176 void ffti( FFT_Tables *fft_tables, double *xr, double *xi, int logm )
  177 {
  178     int nfft = 0;
  179 
  180     kiss_fft_cpx    fin[1 << MAXLOGM];
  181     kiss_fft_cpx    fout[1 << MAXLOGM];
  182 
  183     if ( logm > MAXLOGM )
  184     {
  185         fprintf(stderr, "fft size too big\n");
  186         exit(1);
  187     }
  188 
  189     nfft = logm_to_nfft[logm];
  190 
  191     if ( fft_tables->cfg[logm][1] == NULL )
  192     {
  193         if ( nfft )
  194         {
  195             fft_tables->cfg[logm][1] = kiss_fft_alloc( nfft, 1, NULL, NULL );
  196         }
  197         else
  198         {
  199             fprintf(stderr, "bad logm = %d\n", logm);
  200             exit( 1 );
  201         }
  202     }
  203     
  204     if ( fft_tables->cfg[logm][1] )
  205     {
  206         unsigned int i;
  207         double fac = 1.0 / (double)nfft;
  208         
  209         for ( i = 0; i < nfft; i++ )
  210         {
  211             fin[i].r = xr[i];    
  212             fin[i].i = xi[i];
  213         }
  214         
  215         kiss_fft( (kiss_fft_cfg)fft_tables->cfg[logm][1], fin, fout );
  216 
  217         for ( i = 0; i < nfft; i++ )
  218         {
  219             xr[i]   = fout[i].r * fac;
  220             xi[i]   = fout[i].i * fac;
  221         }
  222     }
  223     else
  224     {
  225         fprintf( stderr, "bad config for logm = %d\n", logm);
  226         exit( 1 );
  227     }
  228 }
  229 
  230 #else /* !defined DRM || defined DRM_1024 */
  231 
  232 void fft_initialize( FFT_Tables *fft_tables )
  233 {
  234     int i;
  235     fft_tables->costbl      = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->costbl[0] ) );
  236     fft_tables->negsintbl   = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->negsintbl[0] ) );
  237     fft_tables->reordertbl  = AllocMemory( (MAXLOGM+1) * sizeof( fft_tables->reordertbl[0] ) );
  238     
  239     for( i = 0; i< MAXLOGM+1; i++ )
  240     {
  241         fft_tables->costbl[i]       = NULL;
  242         fft_tables->negsintbl[i]    = NULL;
  243         fft_tables->reordertbl[i]   = NULL;
  244     }
  245 }
  246 
  247 void fft_terminate( FFT_Tables *fft_tables )
  248 {
  249     int i;
  250 
  251     for( i = 0; i< MAXLOGM+1; i++ )
  252     {
  253         if( fft_tables->costbl[i] != NULL )
  254             FreeMemory( fft_tables->costbl[i] );
  255         
  256         if( fft_tables->negsintbl[i] != NULL )
  257             FreeMemory( fft_tables->negsintbl[i] );
  258             
  259         if( fft_tables->reordertbl[i] != NULL )
  260             FreeMemory( fft_tables->reordertbl[i] );
  261     }
  262 
  263     FreeMemory( fft_tables->costbl );
  264     FreeMemory( fft_tables->negsintbl );
  265     FreeMemory( fft_tables->reordertbl );
  266 
  267     fft_tables->costbl      = NULL;
  268     fft_tables->negsintbl   = NULL;
  269     fft_tables->reordertbl  = NULL;
  270 }
  271 
  272 static void reorder( FFT_Tables *fft_tables, double *x, int logm)
  273 {
  274     int i;
  275     int size = 1 << logm;
  276     unsigned short *r;  //size
  277 
  278 
  279     if ( fft_tables->reordertbl[logm] == NULL ) // create bit reversing table
  280     {
  281         fft_tables->reordertbl[logm] = AllocMemory(size * sizeof(*(fft_tables->reordertbl[0])));
  282 
  283         for (i = 0; i < size; i++)
  284         {
  285             int reversed = 0;
  286             int b0;
  287             int tmp = i;
  288 
  289             for (b0 = 0; b0 < logm; b0++)
  290             {
  291                 reversed = (reversed << 1) | (tmp & 1);
  292                 tmp >>= 1;
  293             }
  294             fft_tables->reordertbl[logm][i] = reversed;
  295         }
  296     }
  297 
  298     r = fft_tables->reordertbl[logm];
  299 
  300     for (i = 0; i < size; i++)
  301     {
  302         int j = r[i];
  303         double tmp;
  304 
  305         if (j <= i)
  306             continue;
  307 
  308         tmp = x[i];
  309         x[i] = x[j];
  310         x[j] = tmp;
  311     }
  312 }
  313 
  314 static void fft_proc(
  315         double *xr, 
  316         double *xi,
  317         fftfloat *refac, 
  318         fftfloat *imfac, 
  319         int size)   
  320 {
  321     int step, shift, pos;
  322     int exp, estep;
  323 
  324     estep = size;
  325     for (step = 1; step < size; step *= 2)
  326     {
  327         int x1;
  328         int x2 = 0;
  329         estep >>= 1;
  330         for (pos = 0; pos < size; pos += (2 * step))
  331         {
  332             x1 = x2;
  333             x2 += step;
  334             exp = 0;
  335             for (shift = 0; shift < step; shift++)
  336             {
  337                 double v2r, v2i;
  338 
  339                 v2r = xr[x2] * refac[exp] - xi[x2] * imfac[exp];
  340                 v2i = xr[x2] * imfac[exp] + xi[x2] * refac[exp];
  341 
  342                 xr[x2] = xr[x1] - v2r;
  343                 xr[x1] += v2r;
  344 
  345                 xi[x2] = xi[x1] - v2i;
  346 
  347                 xi[x1] += v2i;
  348 
  349                 exp += estep;
  350 
  351                 x1++;
  352                 x2++;
  353             }
  354         }
  355     }
  356 }
  357 
  358 static void check_tables( FFT_Tables *fft_tables, int logm)
  359 {
  360     if( fft_tables->costbl[logm] == NULL )
  361     {
  362         int i;
  363         int size = 1 << logm;
  364 
  365         if( fft_tables->negsintbl[logm] != NULL )
  366             FreeMemory( fft_tables->negsintbl[logm] );
  367 
  368         fft_tables->costbl[logm]    = AllocMemory((size / 2) * sizeof(*(fft_tables->costbl[0])));
  369         fft_tables->negsintbl[logm] = AllocMemory((size / 2) * sizeof(*(fft_tables->negsintbl[0])));
  370 
  371         for (i = 0; i < (size >> 1); i++)
  372         {
  373             double theta = 2.0 * M_PI * ((double) i) / (double) size;
  374             fft_tables->costbl[logm][i]     = cos(theta);
  375             fft_tables->negsintbl[logm][i]  = -sin(theta);
  376         }
  377     }
  378 }
  379 
  380 void fft( FFT_Tables *fft_tables, double *xr, double *xi, int logm)
  381 {
  382     if (logm > MAXLOGM)
  383     {
  384         fprintf(stderr, "fft size too big\n");
  385         exit(1);
  386     }
  387 
  388     if (logm < 1)
  389     {
  390         //printf("logm < 1\n");
  391         return;
  392     }
  393 
  394     check_tables( fft_tables, logm);
  395 
  396     reorder( fft_tables, xr, logm);
  397     reorder( fft_tables, xi, logm);
  398 
  399     fft_proc( xr, xi, fft_tables->costbl[logm], fft_tables->negsintbl[logm], 1 << logm );
  400 }
  401 
  402 void rfft( FFT_Tables *fft_tables, double *x, int logm)
  403 {
  404     double xi[1 << MAXLOGR];
  405 
  406     if (logm > MAXLOGR)
  407     {
  408         fprintf(stderr, "rfft size too big\n");
  409         exit(1);
  410     }
  411 
  412     memset(xi, 0, (1 << logm) * sizeof(xi[0]));
  413 
  414     fft( fft_tables, x, xi, logm);
  415 
  416     memcpy(x + (1 << (logm - 1)), xi, (1 << (logm - 1)) * sizeof(*x));
  417 }
  418 
  419 void ffti( FFT_Tables *fft_tables, double *xr, double *xi, int logm)
  420 {
  421     int i, size;
  422     double fac;
  423     double *xrp, *xip;
  424 
  425     fft( fft_tables, xi, xr, logm);
  426 
  427     size = 1 << logm;
  428     fac = 1.0 / size;
  429     xrp = xr;
  430     xip = xi;
  431 
  432     for (i = 0; i < size; i++)
  433     {
  434         *xrp++ *= fac;
  435         *xip++ *= fac;
  436     }
  437 }
  438 
  439 #endif /* defined DRM && !defined DRM_1024 */