real-fft-test.c 4.1 KB

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