-
Notifications
You must be signed in to change notification settings - Fork 1
/
bandpass.go
73 lines (68 loc) · 1.49 KB
/
bandpass.go
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
package biquad
import "math"
type BandPass struct {
blt
}
// NewBandPass creates a band pass filter from
// Fs: sampling frequency
// f0: working frequency
// BW: bandwidth of the filter in octaves
// This filter has constant peak gain (0 dB).
func NewBandPass(Fs, f0, BW float64) (*BandPass, error) {
switch {
case f0 >= Fs:
return nil, ErrBadWorkingFreq
case BW <= 0:
return nil, ErrNegBandwidth
case f0 <= 0 || Fs <= 0:
return nil, ErrBadFreq
}
w0 := 2 * math.Pi * (f0 / Fs)
cos := math.Cos(w0)
alpha := alphaCalc{}.bw(w0, BW)
var (
b0 = alpha
b1 = 0.
b2 = -b0
a0 = 1 + alpha
a1 = -2 * cos
a2 = 1 - alpha
)
return &BandPass{
blt: newBLT(a0, a1, a2, b0, b1, b2),
}, nil
}
// NewBandPass creates a band pass filter from
// Fs: sampling frequency
// Q: peak gain
// BW: bandwidth of the filter in octaves
func NewBandPassFromQ(Fs, Q, BW float64) (*BandPass, error) {
switch {
case BW <= 0:
return nil, ErrNegBandwidth
case Q <= 0:
return nil, ErrBadGain
}
// c = w0/sin(w0).
c := 2 * math.Asinh(1/(2*Q)) / (math.Ln2 * BW)
// apply two newton iterations to obtain w0 starting at sampling frequency.
// f(x) = c*sin(x) - x
// f'(x) = c*cos(x) - 1
w0 := Fs
for i := 0; i < 2; i++ {
w0 -= (c*math.Sin(w0) - w0) / (c*math.Cos(w0) - 1)
}
sin, cos := math.Sincos(w0)
alpha := alphaCalc{}.bw(w0, BW)
var (
b0 = sin / 2
b1 = 0.
b2 = -b0
a0 = 1 + alpha
a1 = -2 * cos
a2 = 1 - alpha
)
return &BandPass{
blt: newBLT(a0, a1, a2, b0, b1, b2),
}, nil
}