A stereo example.
[opus.git] / tests / dft-test.c
1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #include <stdio.h>
6 #include "kiss_fft.h"
7
8 #include "../libcelt/kiss_fft.c"
9
10
11 #ifndef M_PI
12 #define M_PI 3.141592653
13 #endif
14
15 int ret = 0;
16
17 void check(kiss_fft_cpx  * in,kiss_fft_cpx  * out,int nfft,int isinverse)
18 {
19     int bin,k;
20     double errpow=0,sigpow=0, snr;
21     
22     for (bin=0;bin<nfft;++bin) {
23         double ansr = 0;
24         double ansi = 0;
25         double difr;
26         double difi;
27
28         for (k=0;k<nfft;++k) {
29             double phase = -2*M_PI*bin*k/nfft;
30             double re = cos(phase);
31             double im = sin(phase);
32             if (isinverse)
33                 im = -im;
34
35             if (!isinverse)
36             {
37                re /= nfft;
38                im /= nfft;
39             }
40
41             ansr += in[k].r * re - in[k].i * im;
42             ansi += in[k].r * im + in[k].i * re;
43         }
44         /*printf ("%d %d ", (int)ansr, (int)ansi);*/
45         difr = ansr - out[bin].r;
46         difi = ansi - out[bin].i;
47         errpow += difr*difr + difi*difi;
48         sigpow += ansr*ansr+ansi*ansi;
49     }
50     snr = 10*log10(sigpow/errpow);
51     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
52     if (snr<60) {
53        printf( "** poor snr: %f ** \n", snr);
54        ret = 1;
55     }
56 }
57
58 void test1d(int nfft,int isinverse)
59 {
60     size_t buflen = sizeof(kiss_fft_cpx)*nfft;
61
62     kiss_fft_cpx  * in = (kiss_fft_cpx*)malloc(buflen);
63     kiss_fft_cpx  * out= (kiss_fft_cpx*)malloc(buflen);
64     kiss_fft_cfg  cfg = kiss_fft_alloc(nfft,0,0);
65     int k;
66
67     for (k=0;k<nfft;++k) {
68         in[k].r = (rand() % 32767) - 16384;
69         in[k].i = (rand() % 32767) - 16384;
70     }
71
72 #ifdef DOUBLE_PRECISION
73     for (k=0;k<nfft;++k) {
74        in[k].r *= 32768;
75        in[k].i *= 32768;
76     }
77 #endif
78     
79     if (isinverse)
80     {
81        for (k=0;k<nfft;++k) {
82           in[k].r /= nfft;
83           in[k].i /= nfft;
84        }
85     }
86     
87     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
88        
89     if (isinverse)
90        kiss_ifft(cfg,in,out);
91     else
92        kiss_fft(cfg,in,out);
93
94     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
95
96     check(in,out,nfft,isinverse);
97
98     free(in);
99     free(out);
100     free(cfg);
101 }
102
103 int main(int argc,char ** argv)
104 {
105     if (argc>1) {
106         int k;
107         for (k=1;k<argc;++k) {
108             test1d(atoi(argv[k]),0);
109             test1d(atoi(argv[k]),1);
110         }
111     }else{
112         test1d(32,0);
113         test1d(32,1);
114         test1d(128,0);
115         test1d(128,1);
116         test1d(256,0);
117         test1d(256,1);
118 #ifndef RADIX_TWO_ONLY
119         test1d(36,0);
120         test1d(36,1);
121         test1d(50,0);
122         test1d(50,1);
123         test1d(120,0);
124         test1d(120,1);
125         test1d(105,0);
126         test1d(105,1);
127 #endif
128     }
129     return ret;
130 }