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