Changed ordering of real FFT freq data to something more aligned (it's not
[opus.git] / tests / real-fft-test.c
1 #include "kiss_fftr.h"
2 #include "_kiss_fft_guts.h"
3 #include <sys/times.h>
4 #include <time.h>
5 #include <unistd.h>
6 #include <stdio.h>
7 #include <string.h>
8
9 static double cputime(void)
10 {
11     struct tms t;
12     times(&t);
13     return (double)(t.tms_utime + t.tms_stime)/  sysconf(_SC_CLK_TCK) ;
14 }
15
16 static
17 kiss_fft_scalar rand_scalar(void) 
18 {
19 #ifdef USE_SIMD
20     return _mm_set1_ps(rand()-RAND_MAX/2);
21 #else
22     kiss_fft_scalar s = (kiss_fft_scalar)(rand() -RAND_MAX/2);
23     return s/2;
24 #endif
25 }
26
27 static
28 double snr_compare( kiss_fft_cpx * vec1,kiss_fft_scalar * vec2, int n)
29 {
30     int k;
31     double sigpow=1e-10, noisepow=1e-10, err,snr;
32
33     for (k=1;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         noisepow += err * err;
38         err = (double)vec1[k].i - (double)vec2[2*k+1];
39         noisepow += err * err;
40
41     }
42     snr = 10*log10( sigpow / noisepow );
43     if (snr<10) {
44         printf( "\npoor snr: %f\n", snr);
45         exit(1);
46     }
47     return snr;
48 }
49
50 static
51 double snr_compare_scal( kiss_fft_scalar * vec1,kiss_fft_scalar * vec2, int n)
52 {
53     int k;
54     double sigpow=1e-10, noisepow=1e-10, err,snr;
55
56     for (k=1;k<n;++k) {
57         sigpow += (double)vec1[k] * (double)vec1[k];
58         err = (double)vec1[k] - (double)vec2[k];
59         noisepow += err * err;
60     }
61     snr = 10*log10( sigpow / noisepow );
62     if (snr<10) {
63         printf( "\npoor snr: %f\n", snr);
64         exit(1);
65     }
66     return snr;
67 }
68 #define NFFT 8*3*5
69
70 #ifndef NUMFFTS
71 #define NUMFFTS 10000
72 #endif
73
74
75 int main(void)
76 {
77     double ts,tfft,trfft;
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     srand(time(0));
92
93     for (i=0;i<NFFT;++i) {
94         rin[i] = rand_scalar();
95         cin[i].r = rin[i];
96         cin[i].i = zero;
97     }
98
99     kiss_fft_state = kiss_fft_alloc(NFFT,0,0,0);
100     kiss_fftr_state = kiss_fftr_alloc(NFFT,0,0,0);
101     kiss_fft(kiss_fft_state,cin,cout);
102     kiss_fftr(kiss_fftr_state,rin,sout);
103     
104     printf( "nfft=%d, inverse=%d, snr=%g\n",
105             NFFT,0, snr_compare(cout,sout,(NFFT/2)) );
106     ts = cputime();
107     for (i=0;i<NUMFFTS;++i) {
108         kiss_fft(kiss_fft_state,cin,cout);
109     }
110     tfft = cputime() - ts;
111     
112     ts = cputime();
113     for (i=0;i<NUMFFTS;++i) {
114         kiss_fftr( kiss_fftr_state, rin, sout );
115         /* kiss_fftri(kiss_fftr_state,cout,rin); */
116     }
117     trfft = cputime() - ts;
118
119     printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
120
121     free(kiss_fft_state);
122     free(kiss_fftr_state);
123
124     kiss_fft_state = kiss_fft_alloc(NFFT,1,0,0);
125     kiss_fftr_state = kiss_fftr_alloc(NFFT,1,0,0);
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     fin[0] = cin[0].r;
149     fin[1] = cin[NFFT/2].r;
150     for (i=1;i< NFFT/2;++i)
151     {
152        fin[2*i] = cin[i].r;
153        fin[2*i+1] = cin[i].i;
154     }
155     
156     kiss_fft(kiss_fft_state,cin,cout);
157     kiss_fftri(kiss_fftr_state,fin,rout);
158     /*
159     printf(" results from inverse kiss_fft : (%f,%f), (%f,%f), (%f,%f), (%f,%f), (%f,%f) ...\n "
160             , (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
161             ); 
162
163     printf(" results from inverse kiss_fftr: %f,%f,%f,%f,%f ... \n"
164             ,(float)rout[0] ,(float)rout[1] ,(float)rout[2] ,(float)rout[3] ,(float)rout[4]);
165 */
166     for (i=0;i<NFFT;++i) {
167         sout[i] = cout[i].r;
168     }
169
170     printf( "nfft=%d, inverse=%d, snr=%g\n",
171             NFFT,1, snr_compare_scal(rout,sout,NFFT) );
172     free(kiss_fft_state);
173     free(kiss_fftr_state);
174
175     return 0;
176 }