/Bessel function, z-order hyperbolic
zorderHyperbolicBessel(float x)=>
float besselAccuracy = 0.000001
float bessel = 1.0
float summ = 0
float temp = 0
float k = 2.0
float factorial = 1.0
temp := x / 2
summ := temp * temp
bessel += summ
while summ > besselAccuracy
factorial := factorial * k
temp *= x / 2
summ := temp / factorial
summ := summ * summ
bessel += summ
k += 1.0
bessel
//Filter length estimations
filterOrder(float PassBandRipple, float StopBandAttenuation, float PassBandBars, float StopBandBars)=>
float sbripple = 0
float pbripple = 0
float ripple = 0
float attenuation = 0
float bandwidth = 0
float d = 0
float n = 0
float x = 0
float alpha = 0
float FilterLength = 0.
PassBand = 1 / PassBandBars
StopBand = 1 / StopBandBars
bandwidth := PassBand + StopBand
if bandwidth >= 0.5
PassBand := 0.5 * PassBand / bandwidth
StopBand := 0.5 * StopBand / bandwidth
sbripple := math.pow(10.0, (-0.05 * StopBandAttenuation))
pbripple := math.pow(10.0, (0.05 * PassBandRipple)) - 1.0
ripple := math.min(sbripple, pbripple)
attenuation := -20 * math.log(ripple) / math.log(10)
if math.round(attenuation, 5) <= 21.0
alpha := 0.0
d := 0.9222
if math.round(attenuation, 5) > 50.0
alpha := 0.1102 * (attenuation - 8.7)
d := (attenuation - 7.95) / 14.36
if math.round(attenuation, 5) > 21.0 and math.round(attenuation, 5) <= 50
alpha := (0.5842 * math.pow((attenuation - 21.0), 0.4)) + (0.07886 * (attenuation - 21.0))
d := (attenuation - 7.95) / 14.36
n := (d / StopBand) + 1.0
x := math.round(n)
if x % 2 < 1
FilterLength := x
else
FilterLength := x - 1
[int(FilterLength), alpha, PassBand, StopBand]
Normalization(float PassBandRipple, float StopBandAttenuation, float PassBandBars, float StopBandBars)=>
float filter = 0
float Ioalfa = 0
float temp = 0
float norm = 0
[FilterLength, alpha, PassBand, StopBand] = filterOrder(PassBandRipple, StopBandAttenuation, PassBandBars, StopBandBars)
int M = int(FilterLength / 2)
float[] filterCoeff = array.new<float>(FilterLength + 1, 0)
float[] kaiserWindow = array.new<float>(M + 1, 0)
//Window function
norm := M
Ioalfa := zorderHyperbolicBessel(alpha)
for i = 1 to M
temp := i / norm
array.set(kaiserWindow, i, zorderHyperbolicBessel(alpha * math.sqrt(1 - (temp * temp))) / Ioalfa)
//filter coefficients
array.set(filterCoeff, 0, 2.0 * (PassBand + (0.5 * StopBand)))
norm := array.get(filterCoeff, 0)
temp := math.pi * array.get(filterCoeff, 0)
for i = 1 to M
array.set(filterCoeff, i, math.sin(i * temp) * array.get(kaiserWindow, i) / (i * math.pi))
norm := norm + (2 * array.get(filterCoeff, i))
//casual conversion and normalization
float[] NormCoef = array.new<float>(FilterLength + 1, 0)
for i = M + 1 to FilterLength
array.set(filterCoeff, i, array.get(filterCoeff, i - M))
for i = 0 to M - 1
array.set(filterCoeff, i, array.get(filterCoeff, FilterLength - i))
array.set(filterCoeff, M, 2.0 * (PassBand + (0.5 * StopBand)))
for i = 0 to FilterLength
array.set(NormCoef, i, array.get(filterCoeff, i) / norm)
[NormCoef, FilterLength]
filterResponse(float src, float[] NormCoef, int per)=>
float valueBuf = 0
float temp = 0
float temp1 = 0
float Response = 0.0
int i = 0
int filterlength = 0
while filterlength <= per
valueBuf := nz(src[filterlength])
Response := Response + valueBuf * array.get(NormCoef, filterlength)
filterlength += 1
Response