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