Split the radix functions into forward and backward versions, removed the
[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,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     if (isinverse)
54        kiss_ifft(cfg,in,out);
55     else
56        kiss_fft(cfg,in,out);
57
58     check(in,out,nfft,isinverse);
59
60     free(in);
61     free(out);
62     free(cfg);
63 }
64
65 int main(int argc,char ** argv)
66 {
67     if (argc>1) {
68         int k;
69         for (k=1;k<argc;++k) {
70             test1d(atoi(argv[k]),0);
71             test1d(atoi(argv[k]),1);
72         }
73     }else{
74         test1d(32,0);
75         test1d(32,1);
76         test1d(36,0);
77         test1d(36,1);
78         test1d(50,0);
79         test1d(50,1);
80         test1d(120,0);
81         test1d(120,1);
82         test1d(105,0);
83         test1d(105,1);
84     }
85     return 0;
86 }