Split the radix functions into forward and backward versions, removed the
[opus.git] / tests / real-fft-test.c
1 #include "kiss_fftr.h"
2 #include "_kiss_fft_guts.h"
3 #include <sys/times.h>
4 #include <time.h>
5 #include <unistd.h>
6 #include <stdio.h>
7 #include <string.h>
8
9 static double cputime(void)
10 {
11     struct tms t;
12     times(&t);
13     return (double)(t.tms_utime + t.tms_stime)/  sysconf(_SC_CLK_TCK) ;
14 }
15
16 static
17 kiss_fft_scalar rand_scalar(void) 
18 {
19 #ifdef USE_SIMD
20     return _mm_set1_ps(rand()-RAND_MAX/2);
21 #else
22     kiss_fft_scalar s = (kiss_fft_scalar)(rand() -RAND_MAX/2);
23     return s/2;
24 #endif
25 }
26
27 static
28 double snr_compare( kiss_fft_cpx * vec1,kiss_fft_scalar * vec2, int n)
29 {
30     int k;
31     double sigpow=1e-10, noisepow=1e-10, err,snr;
32
33     for (k=1;k<n;++k) {
34         sigpow += (double)vec1[k].r * (double)vec1[k].r + 
35                   (double)vec1[k].i * (double)vec1[k].i;
36         err = (double)vec1[k].r - (double)vec2[2*k];
37         noisepow += err * err;
38         err = (double)vec1[k].i - (double)vec2[2*k+1];
39         noisepow += err * err;
40
41     }
42     snr = 10*log10( sigpow / noisepow );
43     if (snr<10) {
44         printf( "\npoor snr: %f\n", snr);
45         exit(1);
46     }
47     return snr;
48 }
49
50 static
51 double snr_compare_scal( kiss_fft_scalar * vec1,kiss_fft_scalar * vec2, int n)
52 {
53     int k;
54     double sigpow=1e-10, noisepow=1e-10, err,snr;
55
56     for (k=1;k<n;++k) {
57         sigpow += (double)vec1[k] * (double)vec1[k];
58         err = (double)vec1[k] - (double)vec2[k];
59         noisepow += err * err;
60     }
61     snr = 10*log10( sigpow / noisepow );
62     if (snr<10) {
63         printf( "\npoor snr: %f\n", snr);
64         exit(1);
65     }
66     return snr;
67 }
68 #define NFFT 8*3*5
69
70 #ifndef NUMFFTS
71 #define NUMFFTS 10000
72 #endif
73
74
75 int main(void)
76 {
77     double ts,tfft,trfft;
78     int i;
79     kiss_fft_cpx cin[NFFT];
80     kiss_fft_cpx cout[NFFT];
81     kiss_fft_scalar fin[NFFT];
82     kiss_fft_scalar sout[NFFT];
83     kiss_fft_cfg  kiss_fft_state;
84     kiss_fftr_cfg  kiss_fftr_state;
85
86     kiss_fft_scalar rin[NFFT+2];
87     kiss_fft_scalar rout[NFFT+2];
88     kiss_fft_scalar zero;
89     memset(&zero,0,sizeof(zero) ); // ugly way of setting short,int,float,double, or __m128 to zero
90
91     srand(time(0));
92
93     for (i=0;i<NFFT;++i) {
94         rin[i] = rand_scalar();
95         cin[i].r = rin[i];
96         cin[i].i = zero;
97     }
98
99     kiss_fft_state = kiss_fft_alloc(NFFT,0,0);
100     kiss_fftr_state = kiss_fftr_alloc(NFFT,0,0);
101     kiss_fft(kiss_fft_state,cin,cout);
102     kiss_fftr(kiss_fftr_state,rin,sout);
103     
104     printf( "nfft=%d, inverse=%d, snr=%g\n",
105             NFFT,0, snr_compare(cout,sout,(NFFT/2)) );
106     ts = cputime();
107     for (i=0;i<NUMFFTS;++i) {
108         kiss_fft(kiss_fft_state,cin,cout);
109     }
110     tfft = cputime() - ts;
111     
112     ts = cputime();
113     for (i=0;i<NUMFFTS;++i) {
114         kiss_fftr( kiss_fftr_state, rin, sout );
115         /* kiss_fftri(kiss_fftr_state,cout,rin); */
116     }
117     trfft = cputime() - ts;
118
119     printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
120
121     memset(cin,0,sizeof(cin));
122 #if 1
123     cin[0].r = rand_scalar();
124     cin[NFFT/2].r = rand_scalar();
125     for (i=1;i< NFFT/2;++i) {
126         //cin[i].r = (kiss_fft_scalar)(rand()-RAND_MAX/2);
127         cin[i].r = rand_scalar();
128         cin[i].i = rand_scalar();
129     }
130 #else
131     cin[0].r = 12000;
132     cin[3].r = 12000;
133     cin[NFFT/2].r = 12000;
134 #endif
135
136     // conjugate symmetry of real signal 
137     for (i=1;i< NFFT/2;++i) {
138         cin[NFFT-i].r = cin[i].r;
139         cin[NFFT-i].i = - cin[i].i;
140     }
141
142     fin[0] = cin[0].r;
143     fin[1] = cin[NFFT/2].r;
144     for (i=1;i< NFFT/2;++i)
145     {
146        fin[2*i] = cin[i].r;
147        fin[2*i+1] = cin[i].i;
148     }
149     
150     kiss_ifft(kiss_fft_state,cin,cout);
151     kiss_fftri(kiss_fftr_state,fin,rout);
152     /*
153     printf(" results from inverse kiss_fft : (%f,%f), (%f,%f), (%f,%f), (%f,%f), (%f,%f) ...\n "
154             , (float)cout[0].r , (float)cout[0].i , (float)cout[1].r , (float)cout[1].i , (float)cout[2].r , (float)cout[2].i , (float)cout[3].r , (float)cout[3].i , (float)cout[4].r , (float)cout[4].i
155             ); 
156
157     printf(" results from inverse kiss_fftr: %f,%f,%f,%f,%f ... \n"
158             ,(float)rout[0] ,(float)rout[1] ,(float)rout[2] ,(float)rout[3] ,(float)rout[4]);
159 */
160     for (i=0;i<NFFT;++i) {
161         sout[i] = cout[i].r;
162     }
163
164     printf( "nfft=%d, inverse=%d, snr=%g\n",
165             NFFT,1, snr_compare_scal(rout,sout,NFFT) );
166     free(kiss_fft_state);
167     free(kiss_fftr_state);
168
169     return 0;
170 }