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