Made pre-computed twiddles the same for forward and inverse FFT
[opus.git] / tests / dft-test.c
1 #include <stdio.h>
2 #include "kiss_fft.h"
3
4
5 void check(kiss_fft_cpx  * in,kiss_fft_cpx  * out,int nfft,int isinverse)
6 {
7     int bin,k;
8     double errpow=0,sigpow=0;
9     
10     for (bin=0;bin<nfft;++bin) {
11         double ansr = 0;
12         double ansi = 0;
13         double difr;
14         double difi;
15
16         for (k=0;k<nfft;++k) {
17             double phase = -2*M_PI*bin*k/nfft;
18             double re = cos(phase);
19             double im = sin(phase);
20             if (isinverse)
21                 im = -im;
22
23 #ifdef FIXED_POINT
24             re /= nfft;
25             im /= nfft;
26 #endif            
27
28             ansr += in[k].r * re - in[k].i * im;
29             ansi += in[k].r * im + in[k].i * re;
30         }
31         difr = ansr - out[bin].r;
32         difi = ansi - out[bin].i;
33         errpow += difr*difr + difi*difi;
34         sigpow += ansr*ansr+ansi*ansi;
35     }
36     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,10*log10(sigpow/errpow) );
37 }
38
39 void test1d(int nfft,int isinverse)
40 {
41     size_t buflen = sizeof(kiss_fft_cpx)*nfft;
42
43     kiss_fft_cpx  * in = (kiss_fft_cpx*)malloc(buflen);
44     kiss_fft_cpx  * out= (kiss_fft_cpx*)malloc(buflen);
45     kiss_fft_cfg  cfg = kiss_fft_alloc(nfft,isinverse,0,0);
46     int k;
47
48     for (k=0;k<nfft;++k) {
49         in[k].r = (rand() % 65536) - 32768;
50         in[k].i = (rand() % 65536) - 32768;
51     }
52
53     kiss_fft(cfg,in,out);
54
55     check(in,out,nfft,isinverse);
56
57     free(in);
58     free(out);
59     free(cfg);
60 }
61
62 int main(int argc,char ** argv)
63 {
64     if (argc>1) {
65         int k;
66         for (k=1;k<argc;++k) {
67             test1d(atoi(argv[k]),0);
68             test1d(atoi(argv[k]),1);
69         }
70     }else{
71         test1d(32,0);
72         test1d(32,1);
73         test1d(36,0);
74         test1d(36,1);
75         test1d(50,0);
76         test1d(50,1);
77         test1d(120,0);
78         test1d(120,1);
79         test1d(105,0);
80         test1d(105,1);
81     }
82     return 0;
83 }