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