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