lib/celt-0.5.1.3/tests/mdct-test.c

   1 #ifdef HAVE_CONFIG_H
   2 #include "config.h"
   3 #endif
   4
   5 #include <stdio.h>
   6 #include "mdct.c"
   7 #include "stack_alloc.h"
   8 #include "kiss_fft.c"
   9
  10 #ifndef M_PI
  11 #define M_PI 3.141592653
  12 #endif
  13
  14 int ret = 0;
  15 void check(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
  16 {
  17     int bin,k;
  18     double errpow=0,sigpow=0;
  19     double snr;
  20     for (bin=0;bin<nfft/2;++bin) {
  21         double ansr = 0;
  22         double difr;
  23
  24         for (k=0;k<nfft;++k) {
  25            double phase = 2*M_PI*(k+.5+.25*nfft)*(bin+.5)/nfft;
  26            double re = cos(phase);
  27
  28            re /= nfft/4;
  29
  30            ansr += in[k] * re;
  31         }
  32         /*printf ("%f %f\n", ansr, out[bin]);*/
  33         difr = ansr - out[bin];
  34         errpow += difr*difr;
  35         sigpow += ansr*ansr;
  36     }
  37     snr = 10*log10(sigpow/errpow);
  38     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
  39     if (snr<60) {
  40        printf( "** poor snr: %f **\n", snr);
  41        ret = 1;
  42     }
  43 }
  44
  45 void check_inv(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
  46 {
  47    int bin,k;
  48    double errpow=0,sigpow=0;
  49    double snr;
  50    for (bin=0;bin<nfft;++bin) {
  51       double ansr = 0;
  52       double difr;
  53
  54       for (k=0;k<nfft/2;++k) {
  55          double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
  56          double re = cos(phase);
  57
  58          //re *= 2;
  59
  60          ansr += in[k] * re;
  61       }
  62       /*printf ("%f %f\n", ansr, out[bin]);*/
  63       difr = ansr - out[bin];
  64       errpow += difr*difr;
  65       sigpow += ansr*ansr;
  66    }
  67    snr = 10*log10(sigpow/errpow);
  68    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
  69    if (snr<60) {
  70       printf( "** poor snr: %f **\n", snr);
  71       ret = 1;
  72    }
  73 }
  74
  75
  76 void test1d(int nfft,int isinverse)
  77 {
  78     mdct_lookup cfg;
  79     size_t buflen = sizeof(kiss_fft_scalar)*nfft;
  80
  81     kiss_fft_scalar  * in = (kiss_fft_scalar*)malloc(buflen);
  82     kiss_fft_scalar  * out= (kiss_fft_scalar*)malloc(buflen);
  83     celt_word16_t  * window= (celt_word16_t*)malloc(sizeof(celt_word16_t)*nfft/2);
  84     int k;
  85
  86     mdct_init(&cfg, nfft);
  87     for (k=0;k<nfft;++k) {
  88         in[k] = (rand() % 32768) - 16384;
  89     }
  90
  91     for (k=0;k<nfft/2;++k) {
  92        window[k] = Q15ONE;
  93     }
  94 #ifdef DOUBLE_PRECISION
  95     for (k=0;k<nfft;++k) {
  96        in[k] *= 32768;
  97     }
  98 #endif
  99
 100     if (isinverse)
 101     {
 102        for (k=0;k<nfft;++k) {
 103           in[k] /= nfft;
 104        }
 105     }
 106
 107     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
 108
 109     if (isinverse)
 110     {
 111        for (k=0;k<nfft;++k)
 112           out[k] = 0;
 113        mdct_backward(&cfg,in,out, window, nfft/2);
 114        check_inv(in,out,nfft,isinverse);
 115     } else {
 116        mdct_forward(&cfg,in,out,window, nfft/2);
 117        check(in,out,nfft,isinverse);
 118     }
 119     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
 120
 121
 122     free(in);
 123     free(out);
 124     mdct_clear(&cfg);
 125 }
 126
 127 int main(int argc,char ** argv)
 128 {
 129     ALLOC_STACK;
 130     if (argc>1) {
 131         int k;
 132         for (k=1;k<argc;++k) {
 133             test1d(atoi(argv[k]),0);
 134             test1d(atoi(argv[k]),1);
 135         }
 136     }else{
 137         test1d(32,0);
 138         test1d(32,1);
 139         test1d(256,0);
 140         test1d(256,1);
 141         test1d(512,0);
 142         test1d(512,1);
 143 #ifndef RADIX_TWO_ONLY
 144         test1d(40,0);
 145         test1d(40,1);
 146         test1d(56,0);
 147         test1d(56,1);
 148         test1d(120,0);
 149         test1d(120,1);
 150         test1d(240,0);
 151         test1d(240,1);
 152         test1d(480,0);
 153         test1d(480,1);
 154 #endif
 155     }
 156     return ret;
 157 }