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