Changeset 304


Ignore:
Timestamp:
05/02/07 22:55:04 (17 years ago)
Author:
zorg
Message:

Simplify code by using array notation. It should also make the task easier to compilers for auto-vectorization

File:
1 edited

Legend:

Unmodified
Added
Removed
  • libmpc/trunk/libmpcpsy/fft_routines.c

    r137 r304  
    198198PowSpec256 ( const float* x, float* erg )
    199199{
    200     const float*  win = Hann_256;
    201     float*        aix = a;
    202     int           i;
    203 
     200    int i;
    204201    // windowing
    205     i = 256;
    206     while (i--)
    207         *aix++ = *x++ * *win++;
    208 
    209     // perform FFT
    210     rdft ( 256, a, ip, w );
     202    for( i = 0; i < 256; ++i )
     203        a[i] = x[i] * Hann_256[i];
     204
     205    rdft(256, a, ip, w); // perform FFT
    211206
    212207    // calculate power
    213     aix = a;    // reset pointer
    214     i   = 128;
    215     while (i--) {
    216         *erg++ = aix[0]*aix[0] + aix[1]*aix[1];
    217         aix += 2;
    218     }
     208    for( i = 0; i < 128; ++i)
     209        erg[i] = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
    219210}
    220211
     
    224215PowSpec1024 ( const float* x, float* erg )
    225216{
    226     const float*  win = Hann_1024;
    227     float*        aix = a;
    228     int           i;
    229 
    230     i = 1024;                   // windowing
    231     while (i--)
    232         *aix++ = *x++ * *win++;
    233 
    234 //    for (i=0; i<1024; i++)
    235 //        a[i] = Hann_1024[i] * ((i==0 ? 0 : i-512) + 1000);
    236 
    237     rdft ( 1024, a, ip, w );    // perform FFT
    238 
    239     aix = a;                    // calculate power
    240     i   = 512;
    241 
    242 
    243     DECONV;
    244 //    for (i = 0; i <= 512; i++ )
    245 //        printf ("%3u %12.6f %12.6f\n", i, a[i+i], a[i+i+1]);
    246 //    exit(1);
    247     while (i--) {
    248         *erg++ = aix[0]*aix[0] + aix[1]*aix[1];
    249         aix += 2;
    250     }
     217    int i;
     218    // windowing
     219    for( i = 0; i < 1024; ++i )
     220        a[i] = x[i] * Hann_1024[i];
     221
     222    rdft(1024, a, ip, w); // perform FFT
     223
     224    // calculate power
     225    for( i = 0; i < 512; ++i)
     226        erg[i] = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
    251227}
    252228
     
    256232PowSpec2048 ( const float* x, float* erg )
    257233{
    258     const float*  win = Hann_1600;
    259     float*        aix = a;
    260     int           i;
    261 
     234    int i;
    262235    // windowing (only 1600 samples available -> centered in 2048!)
    263     memset ( a     , 0, 224*sizeof(*a) );
    264     aix = a + 224;
    265     i   = 1600;
    266     while (i--)
    267         *aix++ = *x++ * *win++;
    268     memset ( a+1824, 0, 224*sizeof(*a) );
    269 
    270     rdft ( 2048, a, ip, w );    // perform FFT
    271 
    272     aix = a;                    // calculate power
    273     i   = 1024;
    274     while (i--) {
    275         *erg++ = aix[0]*aix[0] + aix[1]*aix[1];
    276         aix += 2;
    277     }
     236    memset(a, 0, 224 * sizeof *a);
     237    for( i = 0; i < 1600; ++i )
     238        a[i+224] = x[i] * Hann_1600[i];
     239    memset(a + 1824, 0, 224 * sizeof *a);
     240
     241    rdft(2048, a, ip, w); // perform FFT
     242
     243    // calculate power
     244    for( i = 0; i < 1024; ++i)
     245        erg[i] = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
    278246}
    279247
     
    283251PolarSpec1024 ( const float* x, float* erg, float* phs )
    284252{
    285     const float*  win = Hann_1024;
    286     float*        aix = a;
    287     int           i;
    288 
    289     i = 1024;                   // windowing
    290     while (i--)
    291         *aix++ = *x++ * *win++;
    292 
    293     rdft ( 1024, a, ip, w );    // perform FFT
     253    int i;
     254    for( i = 0; i < 1024; i++ )
     255        a[i] = x[i] * Hann_1024[i];
     256
     257    rdft( 1024, a, ip, w); // perform FFT
    294258
    295259    // calculate power and phase
    296     aix = a;    // reset pointer
    297     i   = 512;
    298     while (i--) {
    299         *erg++ = aix[0]*aix[0] + aix[1]*aix[1];
    300         *phs++ = ATAN2F (aix[1], aix[0]);
    301         aix += 2;
     260    for( i = 0; i < 512; ++i )
     261    {
     262        erg[i]  = a[i*2] * a[i*2] + a[i*2+1] * a[i*2+1];
     263        phs[i]  = ATAN2F( a[i*2+1], a[i*2] ); 
    302264    }
    303265}
     
    308270Cepstrum2048 ( float* cep, const int MaxLine )
    309271{
    310     float*  aix = cep;
    311     float*  bix = cep + 2048;
    312     int     i;
    313 
     272    int i, j;
    314273    // generate real, even spectrum (symmetric around 1024, cep[2048-i] = cep[i])
    315     for ( i = 0; i < 1024; i++ )
    316         *bix-- = *aix++;
    317 
    318     // perform IFFT
    319     rdft ( 2048, cep, ip, w );
    320 
     274    for( i = 0, j = 1024; i < 1024; ++i, --j )
     275        cep[1024 + j] = cep[i];
     276 
     277    rdft(2048, cep, ip, w);
     278 
    321279    // only real part as outcome (all even indexes of cep[])
    322     aix = cep;
    323     bix = cep;
    324     i   = MaxLine + 1;
    325     while (i--) {
    326         *aix = *bix * (float) (0.9888 / 2048.);
    327 //      *aix = *bix * 0.0004828125f;
    328         aix ++;
    329         bix += 2;
    330     }
    331 }
     280    for( i = 0; i < MaxLine + 1; ++i )
     281        cep[i] = cep[i*2] * (float) (0.9888 / 2048);
     282}
Note: See TracChangeset for help on using the changeset viewer.