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