101c1e835b2dd4a5f11a74766476c95c5760f75d
[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     if (isinverse)
62     {
63        for (k=0;k<nfft;++k) {
64           in[k].r /= nfft;
65           in[k].i /= nfft;
66        }
67     }
68     
69     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
70        
71     if (isinverse)
72        kiss_ifft(cfg,in,out);
73     else
74        kiss_fft(cfg,in,out);
75
76     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
77
78     check(in,out,nfft,isinverse);
79
80     free(in);
81     free(out);
82     free(cfg);
83 }
84
85 int main(int argc,char ** argv)
86 {
87     if (argc>1) {
88         int k;
89         for (k=1;k<argc;++k) {
90             test1d(atoi(argv[k]),0);
91             test1d(atoi(argv[k]),1);
92         }
93     }else{
94         test1d(32,0);
95         test1d(32,1);
96         test1d(36,0);
97         test1d(36,1);
98         test1d(50,0);
99         test1d(50,1);
100         test1d(120,0);
101         test1d(120,1);
102         test1d(105,0);
103         test1d(105,1);
104     }
105     return 0;
106 }