1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
|
#ifndef __FIRfilter
#include "FIRfilter.h"
#endif
//-----------------------------------------------------------------------------
// you're supposed to use use an odd number of taps
void calculateFIRidealLowpassCoefficients(float cutoff, float samplerate,
int numTaps, float *coefficients)
{
int middleCoeff;
float corner, value;
// get the cutoff as a ratio of cutoff to Nyquist, scaled from 0 to Pi
corner = (cutoff / (samplerate*0.5f)) * PI;
if (numTaps%2)
{
middleCoeff = (numTaps-1) / 2;
coefficients[middleCoeff] = corner/PI;
}
else
middleCoeff = numTaps/2;
for (int n=0; n < middleCoeff; n++)
{
value = (float)n - ((float)(numTaps-1) * 0.5f);
coefficients[n] = sinf(value*corner) / (value*PI);
coefficients[numTaps-1-n] = coefficients[n];
}
}
//-----------------------------------------------------------------------------
void applyKaiserWindow(int numTaps, float *coefficients, float attenuation)
{
int halfLength;
// beta is 0 if the attenuation is less than 21 dB
float beta = 0.0f;
if (attenuation >= 50.0f)
beta = 0.1102f * (attenuation - 8.71f);
else if ( (attenuation < 50.0f) && (attenuation >= 21.0f) )
{
beta = 0.5842f * powf( (attenuation - 21.0f), 0.4f);
beta += 0.07886f * (attenuation - 21.0f);
}
if (numTaps%2)
halfLength = (numTaps+1) / 2;
else
halfLength = numTaps / 2;
for (int n=0; n < halfLength; n++)
{
coefficients[n] *= besselIzero(beta *
sqrtf(1.0f - powf( (1.0f-((2.0f*n)/((float)(numTaps-1)))), 2.0f )))
/ besselIzero(beta);
coefficients[numTaps-1-n] = coefficients[n];
}
}
//-----------------------------------------------------------------------------
float besselIzero(float in)
{
float sum, numerator, denominator, term, halfIn;
sum = 1.0f;
halfIn = in * 0.5f;
denominator = 1.0f;
numerator = 1.0f;
for (int m=1; m <= 32; m++)
{
numerator *= halfIn;
denominator *= (float)m;
term = numerator / denominator;
sum += term * term;
}
return sum;
}
//-----------------------------------------------------------------------------
float besselIzero2(float in)
{
float sum = 1.0f;
float ds = 1.0f;
float d = 0.0f;
do
{
d += 2.0f;
ds *= ( in * in ) / ( d * d );
sum += ds;
}
while ( ds > (1E-7f * sum) );
return sum;
}
|
