/project/jni/sdl_mixer/timidity/filter.c

https://github.com/aichunyu/FFPlayer · C · 203 lines · 117 code · 30 blank · 56 comment · 20 complexity · 7da1f19b11388bcbd68f528425252bf7 MD5 · raw file

  1. /*
  2. TiMidity -- Experimental MIDI to WAVE converter
  3. Copyright (C) 1995 Tuukka Toivonen <toivonen@clinet.fi>
  4. This program is free software; you can redistribute it and/or modify
  5. it under the terms of the GNU General Public License as published by
  6. the Free Software Foundation; either version 2 of the License, or
  7. (at your option) any later version.
  8. This program is distributed in the hope that it will be useful,
  9. but WITHOUT ANY WARRANTY; without even the implied warranty of
  10. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  11. GNU General Public License for more details.
  12. You should have received a copy of the GNU General Public License
  13. along with this program; if not, write to the Free Software
  14. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  15. filter.c: written by Vincent Pagel ( pagel@loria.fr )
  16. implements fir antialiasing filter : should help when setting sample
  17. rates as low as 8Khz.
  18. April 95
  19. - first draft
  20. 22/5/95
  21. - modify "filter" so that it simulate leading and trailing 0 in the buffer
  22. */
  23. #include <stdio.h>
  24. #include <string.h>
  25. #include <math.h>
  26. #include <stdlib.h>
  27. #include "config.h"
  28. #include "common.h"
  29. #include "ctrlmode.h"
  30. #include "instrum.h"
  31. #include "filter.h"
  32. /* bessel function */
  33. static float ino(float x)
  34. {
  35. float y, de, e, sde;
  36. int i;
  37. y = x / 2;
  38. e = 1.0;
  39. de = 1.0;
  40. i = 1;
  41. do {
  42. de = de * y / (float) i;
  43. sde = de * de;
  44. e += sde;
  45. } while (!( (e * 1.0e-08 - sde > 0) || (i++ > 25) ));
  46. return(e);
  47. }
  48. /* Kaiser Window (symetric) */
  49. static void kaiser(float *w,int n,float beta)
  50. {
  51. float xind, xi;
  52. int i;
  53. xind = (float)((2*n - 1) * (2*n - 1));
  54. for (i =0; i<n ; i++)
  55. {
  56. xi = (float)(i + 0.5);
  57. w[i] = ino((float)(beta * sqrt((double)(1. - 4 * xi * xi / xind))))
  58. / ino((float)beta);
  59. }
  60. }
  61. /*
  62. * fir coef in g, cuttoff frequency in fc
  63. */
  64. static void designfir(float *g , float fc)
  65. {
  66. int i;
  67. float xi, omega, att, beta ;
  68. float w[ORDER2];
  69. for (i =0; i < ORDER2 ;i++)
  70. {
  71. xi = (float) (i + 0.5);
  72. omega = (float)(PI * xi);
  73. g[i] = (float)(sin( (double) omega * fc) / omega);
  74. }
  75. att = 40.; /* attenuation in db */
  76. beta = (float) (exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886
  77. * (att - 20.96));
  78. kaiser( w, ORDER2, beta);
  79. /* Matrix product */
  80. for (i =0; i < ORDER2 ; i++)
  81. g[i] = g[i] * w[i];
  82. }
  83. /*
  84. * FIR filtering -> apply the filter given by coef[] to the data buffer
  85. * Note that we simulate leading and trailing 0 at the border of the
  86. * data buffer
  87. */
  88. static void filter(sample_t *result,sample_t *data, int32 length,float coef[])
  89. {
  90. int32 sample,i,sample_window;
  91. int16 peak = 0;
  92. float sum;
  93. /* Simulate leading 0 at the begining of the buffer */
  94. for (sample = 0; sample < ORDER2 ; sample++ )
  95. {
  96. sum = 0.0;
  97. sample_window= sample - ORDER2;
  98. for (i = 0; i < ORDER ;i++)
  99. sum += (float)(coef[i] *
  100. ((sample_window<0)? 0.0 : data[sample_window++])) ;
  101. /* Saturation ??? */
  102. if (sum> 32767.) { sum=32767.; peak++; }
  103. if (sum< -32768.) { sum=-32768; peak++; }
  104. result[sample] = (sample_t) sum;
  105. }
  106. /* The core of the buffer */
  107. for (sample = ORDER2; sample < length - ORDER + ORDER2 ; sample++ )
  108. {
  109. sum = 0.0;
  110. sample_window= sample - ORDER2;
  111. for (i = 0; i < ORDER ;i++)
  112. sum += data[sample_window++] * coef[i];
  113. /* Saturation ??? */
  114. if (sum> 32767.) { sum=32767.; peak++; }
  115. if (sum< -32768.) { sum=-32768; peak++; }
  116. result[sample] = (sample_t) sum;
  117. }
  118. /* Simulate 0 at the end of the buffer */
  119. for (sample = length - ORDER + ORDER2; sample < length ; sample++ )
  120. {
  121. sum = 0.0;
  122. sample_window= sample - ORDER2;
  123. for (i = 0; i < ORDER ;i++)
  124. sum += (float)(coef[i] *
  125. ((sample_window>=length)? 0.0 : data[sample_window++])) ;
  126. /* Saturation ??? */
  127. if (sum> 32767.) { sum=32767.; peak++; }
  128. if (sum< -32768.) { sum=-32768; peak++; }
  129. result[sample] = (sample_t) sum;
  130. }
  131. if (peak)
  132. ctl->cmsg(CMSG_ERROR, VERB_NORMAL,
  133. "Saturation %2.3f %%.", 100.0*peak/ (float) length);
  134. }
  135. /***********************************************************************/
  136. /* Prevent aliasing by filtering any freq above the output_rate */
  137. /* */
  138. /* I don't worry about looping point -> they will remain soft if they */
  139. /* were already */
  140. /***********************************************************************/
  141. void antialiasing(Sample *sp, int32 output_rate )
  142. {
  143. sample_t *temp;
  144. int i;
  145. float fir_symetric[ORDER];
  146. float fir_coef[ORDER2];
  147. float freq_cut; /* cutoff frequency [0..1.0] FREQ_CUT/SAMP_FREQ*/
  148. ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: Fsample=%iKHz",
  149. sp->sample_rate);
  150. /* No oversampling */
  151. if (output_rate>=sp->sample_rate)
  152. return;
  153. freq_cut= (float) output_rate / (float) sp->sample_rate;
  154. ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: cutoff=%f%%",
  155. freq_cut*100.);
  156. designfir(fir_coef,freq_cut);
  157. /* Make the filter symetric */
  158. for (i = 0 ; i<ORDER2 ;i++)
  159. fir_symetric[ORDER-1 - i] = fir_symetric[i] = fir_coef[ORDER2-1 - i];
  160. /* We apply the filter we have designed on a copy of the patch */
  161. temp = safe_malloc(sp->data_length);
  162. memcpy(temp,sp->data,sp->data_length);
  163. filter(sp->data,temp,sp->data_length/sizeof(sample_t),fir_symetric);
  164. free(temp);
  165. }