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