kaiser.c
上传用户:dangjiwu
上传日期:2013-07-19
资源大小:42019k
文件大小:5k
源码类别:

Symbian

开发平台:

Visual C++

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