Mode doc
[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 <sys/times.h>
8 #include <time.h>
9 #include <unistd.h>
10 #include <stdio.h>
11 #include <string.h>
12
13 static double cputime(void)
14 {
15     struct tms t;
16     times(&t);
17     return (double)(t.tms_utime + t.tms_stime)/  sysconf(_SC_CLK_TCK) ;
18 }
19
20 static
21 kiss_fft_scalar rand_scalar(void) 
22 {
23 #ifdef USE_SIMD
24     return _mm_set1_ps(rand()-RAND_MAX/2);
25 #else
26     kiss_fft_scalar s = (kiss_fft_scalar)(rand() -RAND_MAX/2);
27     return s/2;
28 #endif
29 }
30
31 static
32 double snr_compare( kiss_fft_cpx * vec1,kiss_fft_scalar * vec2, int n)
33 {
34     int k;
35     double sigpow=1e-10, noisepow=1e-10, err,snr;
36
37     for (k=1;k<n;++k) {
38         sigpow += (double)vec1[k].r * (double)vec1[k].r + 
39                   (double)vec1[k].i * (double)vec1[k].i;
40         err = (double)vec1[k].r - (double)vec2[2*k];
41         noisepow += err * err;
42         err = (double)vec1[k].i - (double)vec2[2*k+1];
43         noisepow += err * err;
44
45     }
46     snr = 10*log10( sigpow / noisepow );
47     if (snr<10) {
48         printf( "\npoor snr: %f\n", snr);
49         exit(1);
50     }
51     return snr;
52 }
53
54 static
55 double snr_compare_scal( kiss_fft_scalar * vec1,kiss_fft_scalar * vec2, int n)
56 {
57     int k;
58     double sigpow=1e-10, noisepow=1e-10, err,snr;
59
60     for (k=1;k<n;++k) {
61         sigpow += (double)vec1[k] * (double)vec1[k];
62         err = (double)vec1[k] - (double)vec2[k];
63         noisepow += err * err;
64     }
65     snr = 10*log10( sigpow / noisepow );
66     if (snr<10) {
67         printf( "\npoor snr: %f\n", snr);
68         exit(1);
69     }
70     return snr;
71 }
72 #define NFFT 8*3*5
73
74 #ifndef NUMFFTS
75 #define NUMFFTS 10000
76 #endif
77
78
79 int main(void)
80 {
81     double ts,tfft,trfft;
82     int i;
83     kiss_fft_cpx cin[NFFT];
84     kiss_fft_cpx cout[NFFT];
85     kiss_fft_scalar fin[NFFT];
86     kiss_fft_scalar sout[NFFT];
87     kiss_fft_cfg  kiss_fft_state;
88     kiss_fftr_cfg  kiss_fftr_state;
89
90     kiss_fft_scalar rin[NFFT+2];
91     kiss_fft_scalar rout[NFFT+2];
92     kiss_fft_scalar zero;
93     memset(&zero,0,sizeof(zero) ); // ugly way of setting short,int,float,double, or __m128 to zero
94
95     srand(time(0));
96
97     for (i=0;i<NFFT;++i) {
98         rin[i] = rand_scalar();
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     ts = cputime();
111     for (i=0;i<NUMFFTS;++i) {
112         kiss_fft(kiss_fft_state,cin,cout);
113     }
114     tfft = cputime() - ts;
115     
116     ts = cputime();
117     for (i=0;i<NUMFFTS;++i) {
118         kiss_fftr( kiss_fftr_state, rin, sout );
119         /* kiss_fftri(kiss_fftr_state,cout,rin); */
120     }
121     trfft = cputime() - ts;
122
123     printf("%d complex ffts took %gs, real took %gs\n",NUMFFTS,tfft,trfft);
124
125     memset(cin,0,sizeof(cin));
126 #if 1
127     cin[0].r = rand_scalar();
128     cin[NFFT/2].r = rand_scalar();
129     for (i=1;i< NFFT/2;++i) {
130         //cin[i].r = (kiss_fft_scalar)(rand()-RAND_MAX/2);
131         cin[i].r = rand_scalar();
132         cin[i].i = rand_scalar();
133     }
134 #else
135     cin[0].r = 12000;
136     cin[3].r = 12000;
137     cin[NFFT/2].r = 12000;
138 #endif
139
140     // conjugate symmetry of real signal 
141     for (i=1;i< NFFT/2;++i) {
142         cin[NFFT-i].r = cin[i].r;
143         cin[NFFT-i].i = - cin[i].i;
144     }
145
146     fin[0] = cin[0].r;
147     fin[1] = cin[NFFT/2].r;
148     for (i=1;i< NFFT/2;++i)
149     {
150        fin[2*i] = cin[i].r;
151        fin[2*i+1] = cin[i].i;
152     }
153     
154     kiss_ifft(kiss_fft_state,cin,cout);
155     kiss_fftri(kiss_fftr_state,fin,rout);
156     /*
157     printf(" results from inverse kiss_fft : (%f,%f), (%f,%f), (%f,%f), (%f,%f), (%f,%f) ...\n "
158             , (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
159             ); 
160
161     printf(" results from inverse kiss_fftr: %f,%f,%f,%f,%f ... \n"
162             ,(float)rout[0] ,(float)rout[1] ,(float)rout[2] ,(float)rout[3] ,(float)rout[4]);
163 */
164     for (i=0;i<NFFT;++i) {
165         sout[i] = cout[i].r;
166     }
167
168     printf( "nfft=%d, inverse=%d, snr=%g\n",
169             NFFT,1, snr_compare_scal(rout,sout,NFFT) );
170     free(kiss_fft_state);
171     free(kiss_fftr_state);
172
173     return 0;
174 }