kaiser.c
上传用户:zhongxx05
上传日期:2007-06-06
资源大小:33641k
文件大小:4k
源码类别:

Symbian

开发平台:

C/C++

  1. /* ***** BEGIN LICENSE BLOCK ***** 
  2.  * Version: RCSL 1.0/RPSL 1.0 
  3.  *  
  4.  * Portions Copyright (c) 1995-2002 RealNetworks, Inc. All Rights Reserved. 
  5.  *      
  6.  * The contents of this file, and the files included with this file, are 
  7.  * subject to the current version of the RealNetworks Public Source License 
  8.  * Version 1.0 (the "RPSL") available at 
  9.  * http://www.helixcommunity.org/content/rpsl unless you have licensed 
  10.  * the file under the RealNetworks Community Source License Version 1.0 
  11.  * (the "RCSL") available at http://www.helixcommunity.org/content/rcsl, 
  12.  * in which case the RCSL will apply. You may also obtain the license terms 
  13.  * directly from RealNetworks.  You may not use this file except in 
  14.  * compliance with the RPSL or, if you have a valid RCSL with RealNetworks 
  15.  * applicable to this file, the RCSL.  Please see the applicable RPSL or 
  16.  * RCSL for the rights, obligations and limitations governing use of the 
  17.  * contents of the file.  
  18.  *  
  19.  * This file is part of the Helix DNA Technology. RealNetworks is the 
  20.  * developer of the Original Code and owns the copyrights in the portions 
  21.  * it created. 
  22.  *  
  23.  * This file, and the files included with this file, is distributed and made 
  24.  * available on an 'AS IS' basis, WITHOUT WARRANTY OF ANY KIND, EITHER 
  25.  * EXPRESS OR IMPLIED, AND REALNETWORKS HEREBY DISCLAIMS ALL SUCH WARRANTIES, 
  26.  * INCLUDING WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY, FITNESS 
  27.  * FOR A PARTICULAR PURPOSE, QUIET ENJOYMENT OR NON-INFRINGEMENT. 
  28.  * 
  29.  * Technology Compatibility Kit Test Suite(s) Location: 
  30.  *    http://www.helixcommunity.org/content/tck 
  31.  * 
  32.  * Contributor(s): 
  33.  *  
  34.  * ***** END LICENSE BLOCK ***** */ 
  35. #include "hlxclib/math.h"
  36. #include "kaiser.h"
  37. #include "allresamplers.h"
  38. #ifndef M_PI
  39. #define M_PI 3.14159265358979323846
  40. #endif
  41. #define IZEROEPS (1E-21) /* max error in Izero */
  42. /*
  43.  * KaiserEstim() estimates the window length and beta needed to
  44.  * meet the given filter specifications.
  45.  *
  46.  * fpass, fstop are normalized freq (1.0 == Nyquist).
  47.  * atten is stopband attenuation in dB.
  48.  */
  49. void
  50. KaiserEstim(float fpass, float fstop, float atten, int *length, float *beta)
  51. {
  52. double d, b;
  53. /* estimate the required beta */
  54. if (atten < 21.0f) {
  55. b = 0.0;
  56. } else if (atten <= 50.0f) {
  57. b = 0.5842 * pow((atten-21.0), 0.4) + 0.07886 * (atten-21.0);
  58. } else { /* atten > 50 */
  59. b = 0.1102 * (atten-8.7);
  60. }
  61. *beta = (float)b;
  62. /* estimate the required length */
  63. d = (atten - 7.95) / (M_PI * 2.285);
  64. *length = 1 + (int)(d / (fstop - fpass));
  65. }
  66. /*
  67.  *             inf
  68.  * Io(x) = 1 + sum(((x/2)^r / r!)^2)
  69.  *             r=1
  70.  */
  71. double
  72. Izero(double x)
  73. {
  74. double halfx, term, sum, r;
  75. double temp;
  76. halfx = 0.5 * x;
  77. term = sum = r = 1.0;
  78. do {
  79. temp = halfx / r;
  80. term *= temp * temp;
  81. sum += term;
  82. r += 1.0;
  83. } while (term >= (sum * IZEROEPS));
  84. return sum;
  85. }
  86. /*
  87.  * KaiserLowpass() creates a Kaiser-windowed lowpass filter.
  88.  *
  89.  * length is length of filter wing.
  90.  * cutoff is normalized cutoff freq (-6dB down).
  91.  * beta is the Kaiser window parameter.
  92.  * gain is the desired dc gain.
  93.  *
  94.  * The Kaiser window is given by:
  95.  * w[n] = Io(beta * sqrt(1 - (n/M)^2)) / Io(beta) (0 <= n <= M)
  96.  *
  97.  * Note: sampling of the "analog" lowpass is offset by 0.5 sample,
  98.  * so that all phases are an even length.
  99.  */
  100. void 
  101. KaiserLowpass(int length, float cutoff, float beta, float gain, double *filter)
  102. {
  103. double ibeta, ilength, x, w;
  104. int i;
  105. ibeta = 1.0 / Izero(beta);
  106. ilength = 1.0 / (length - 0.5); /* last window value should be ibeta */
  107. for (i = 0; i < length; i++) {
  108. x = i + 0.5;
  109. /* Kaiser window */
  110. w = x * ilength;
  111. w = 1.0 - w * w;
  112. w = MAX(w, 0.0);
  113. w = Izero(beta * sqrt(w)) * ibeta;
  114. /* windowed ideal lowpass filter */
  115. filter[i] = w * gain * sin(cutoff * M_PI * x) / (M_PI * x);
  116. }
  117. }