|
 Do a World Wide Web (WWW) search for filter equations
and you will find thousands of pages, including a few here on RF Cafe. However, if you
want an example of how to implement the transfer functions in a spreadsheet or software,
examples of actual code are elusive (other than maybe a Matlab or MathCAD worksheet).
As one who has incorporated equations for Butterworth, Chebyshev Type 1, Chebyshev
Type 2, Bessel, and other filter functions in many spreadsheets and software over
the past few decades, I figured it might be useful to post snippets of my code so that
someone else can copy and paste it directly into other work. BTW, I do not consider myself
to be a filter expert by any means and there is no ground-breaking knowledge here; it's
just hopefully easier to find. Please look at my Filter Transfer
Functions page for some of the textbook equations.
Writing a macro to use in a spreadsheet is the preferred approach for the sake of
maintenance. Excel still supports Visual Basic for Applications (VBA), and since its
code format is very readable, that is what is shown here. They can be adapted very easily
for use in other languages. Since the VBA code below was copied and pasted directly from
my Excel (2007) spreadsheet, you should be able to copy and paste it directly into yours.
The following definitions are used in all the equations:
f = frequency for calculation of gain
| fL = lower cutoff frequency
| fU = upper cutoff frequency
fRel = frequency relative to the cutoff frequency |
BW = bandwidth | N = filter order*
Ripple = gain ripple (dB) |
PBRipple = passband ripple (dB) |
SBRipple = stopband ripple (dB)
* Note that for calculation purposes N does not have to be an integer. That
allows you to tweak the filter model to match the actual measured results. It will have
a noticeable effect on the appearance of inband and out-of-band ripple spacing, but not
min/max amplitude.
 Macro formula calls from within a spreadsheet
cell look similar to this:
=Cheby1BPGain($B48, $C$8, $C$7, $C$16,$C$18) ,
where the Column/Row cell references are where the defined values of:
cell B48 = Frequency for which gain is being calculated
cell C8 = Lower cutoff frequency
cell C7 = Upper cutoff frequency
cell C16 = Filter order
cell C18 = Passband ripple (dB)
Filter function names are colored red in the declaration
and everywhere they are assigned a value.
Since logarithmic charts are used, a check is made to assure negative and zero (0)
values for frequency are not calculated, which could result in a calculation and/or graph
scaling error.
Butterworth Filters | Chebyshev
Type 1 Filters | Chebyshev Type 2 Filters |
Bessel Filters
Butterworth filters are probably the simplest to implement because the equations are
the same both inside and outside the passband.
Public Function ButterLowpassGain(f As
Double, fU As Double, N As Double) As Double
Dim fRel As Double
fRel = f / fU
ButterLowpassGain = -10 * WorksheetFunction.Log10(1
+ (fRel) ^ (2 * N))
End Function
Public Function ButterHighpassGain(f
As Double, fL As Double, N As Double) As Double
Dim fRel As Double
fRel = fL / f
ButterHighpassGain = -10 * WorksheetFunction.Log10(1
+ (fRel) ^ (2 * N))
End Function
Public Function ButterBandpassGain(f
As Double, fL As Double, fU As Double, N As Double) As Double
Dim fRel As Double
Dim fo As Double
Dim BW As Double
fo = (fU * fL) ^ 0.5
BW = (fU - fL) / fo
fRel = ((f / fo) - (fo / f)) / BW
ButterBandpassGain = -10 * WorksheetFunction.Log10(1
+ (fRel) ^ (2 * N))
End Function
Public Function ButterBandstopGain(f
As Double, fL As Double, fU As Double, N As Double) As Double
Dim fRel As Double
Dim fo As Double
Dim BW As Double
fo = (fU * fL) ^ 0.5
BW = (fU - fL) / fo
fRel = 1 / (((f / fo) - (fo / f)) / BW)
ButterBandstopGain = -10 * WorksheetFunction.Log10(1
+ (fRel) ^ (2 * N))
End Function
Chebyshev Type 1 filters have two distinct regions where the transfer function are
different. The inband region is a standard cosine function whereas the out-of-band region
is a hyperbolic cosine function. This requires checking to determine whether the frequency
used for calculation is in-band or out-of-band. Note that the "esqr" term is different
than that of the Chebyshev Type 2 filter.
Public Function Cheby1LPGain(f
As Double, fU As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
fRel = f / fU
esqr = 10 ^ (Ripple / 10) - 1
If fRel < 1 Then 'In-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
Else 'Out-of-band
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
End If
Cheby1LPGain = -10 * WorksheetFunction.Log10(1 +
esqr * Cn ^ 2)
End Function
Public Function Cheby1HPGain(f
As Double, fL As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
fRel = fL / f
esqr = 10 ^ (Ripple / 10) - 1
If fRel < 1 Then 'Out-of-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
Else 'In-band
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
End If
Cheby1HPGain = -10 * WorksheetFunction.Log10(1 +
esqr * Cn ^ 2)
End Function
Public Function Cheby1BPGain(f
As Double, fL As Double, fU As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
Dim fo As Double
Dim BW As Double
fo = Sqr(fU * fL)
BW = (fU - fL) / fo
fRel = Abs(((f / fo) - (fo / f)) / BW)
esqr = 10 ^ (Ripple / 10) - 1
If (f < fL) Or (f > fU) Then
'Out-of-band
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
Else 'In-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
End If
Cheby1BPGain = -10 * WorksheetFunction.Log10(1 +
esqr * Cn ^ 2)
End Function
Public Function Cheby1BSGain(f
As Double, fL As Double, fU As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
Dim fo As Double
Dim BW As Double
fo = Sqr(fU * fL)
BW = (fU - fL) / fo
fRel = 1 / Abs(((f / fo) - (fo / f)) / BW)
esqr = 10 ^ (Ripple / 10) - 1
If (f <= fL) Or (f >= fU) Then
'In-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
Else 'Out-of-band
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
End If
Cheby1BSGain = -10 * WorksheetFunction.Log10(1 +
esqr * Cn ^ 2)
End Function
Chebyshev Type 2 filters have two distinct regions where the transfer function are
different. The inband region is a standard cosine function whereas the out-of-band region
is a hyperbolic cosine function. This requires checking to determine whether the frequency
used for calculation is in-band or out-of-band. Note that the "esqr" term is different
than that of the Chebyshev Type 1 filter.
Public Function Cheby2LPGain(f
As Double, fU As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
fRel = fU / f
esqr = -10 ^ (Ripple / 10)
esqr = esqr / (1 - esqr ^ 2)
If fRel < 1 Then 'In-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
Else 'Out-of-band
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
End If
Cheby2LPGain = -10 * WorksheetFunction.Log10((1 +
esqr * Cn ^ 2) / (esqr * Cn ^ 2))
End Function
Public Function Cheby2HPGain(f
As Double, fL As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
fRel = f / fL
esqr = -10 ^ (Ripple / 10)
esqr = esqr / (1 - esqr ^ 2)
If fRel < 1 Then ' Out-of-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
Else ' Inband
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
End If
Cheby2HPGain = -10 * WorksheetFunction.Log10((1 +
esqr * Cn ^ 2) / (esqr * Cn ^ 2))
End Function
Public Function Cheby2BPGain(f
As Double, fL As Double, fU As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
Dim fo As Double
Dim BW As Double
fo = Sqr(fU * fL)
BW = (fU - fL) / fo
fRel = Abs(((f / fo) - (fo / f)) / BW)
fRel = 1 / fRel
esqr = -10 ^ (Ripple / 10)
esqr = esqr / (1 - esqr ^ 2)
If (f < fL) Or (f > fU) Then
'Out-of-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
Else 'In-band
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
End If
Cheby2BPGain = -10 * WorksheetFunction.Log10((1 +
esqr * Cn ^ 2) / (esqr * Cn ^ 2))
End Function
Public Function Cheby2BSGain(f
As Double, fL As Double, fU As Double, N As Double, Ripple As Double) As Double
Dim fRel As Double
Dim esqr As Double
Dim Cn As Double
Dim fo As Double
Dim BW As Double
fo = Sqr(fU * fL)
BW = (fU - fL) / fo
fRel = Abs(((f / fo) - (fo / f)) / BW)
esqr = -10 ^ (Ripple / 10)
esqr = esqr / (1 - esqr ^ 2)
If (f < fL) Or (f > fU) Then
'Out-of-band
Cn = WorksheetFunction.Cosh(N * WorksheetFunction.Acosh(fRel))
Else 'In-band
Cn = Cos(N * WorksheetFunction.Acos(fRel))
End If
Cheby2BSGain = -10 * WorksheetFunction.Log10((1 +
esqr * Cn ^ 2) / (esqr * Cn ^ 2))
End Function
Bessel filters exhibit lower out-of-band attenuation rates than Butterworth and Chebyshev,
but they have very flat inband linear phase and flat group delay. Textbook suggest a
general normalizing factor of √(2N - 1) * ln2) to place
the 3 dB corner frequency at the specified cutoff frequency. My experimentation
has found that for filter orders 3 through 15, √(1.281.3N - 1) comes closer
to 3 dB at the cutoff frequency. Use whatever factor you are most comfortable with.
Public Function BesselLPGain(f As Double,
fU As Double, N As Double) As Double
Dim Im(0 To 25) As Double
Dim Re(0 To 25) As Double
Dim Bn As Double
Dim Pn As Double
Dim c(0 To 25) As Double
Dim fRel As Double
Dim fAdj As Double
Dim k As Integer
fAdj = (1.28 ^ 1.3 * N - 1) ^ 0.5
fRel = f / fU * fAdj
Bn = 0
Pn = 0
For k = 0 To N
c(k) = WorksheetFunction.Fact(2 * N - k) / _
(2 ^ (N - k) *
WorksheetFunction.Fact(k) * WorksheetFunction.Fact(N - k))
If WorksheetFunction.IsOdd(k) Then
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(Abs(k -
1) / 2, 0)) Then
Im(k) = -c(k)
Else
Im(k) = c(k)
End If
Pn = Pn + fRel ^ k * Im(k)
Else
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(k / 2, 0))
Then
Re(k) = -c(k)
Else
Re(k) = c(k)
End If
Bn = Bn + fRel ^ k * Re(k)
End If
Next k
BesselLPGain = 20 * WorksheetFunction.Log10(c(0)
/ (Bn ^ 2 + Pn ^ 2) ^ 0.5)
End Function
Public Function BesselHPGain(f As Double,
fL As Double, N As Double) As Double
Dim Im(0 To 25) As Double
Dim Re(0 To 25) As Double
Dim Bn As Double
Dim Pn As Double
Dim c(0 To 25) As Double
Dim fRel As Double
Dim fAdj As Double
Dim k As Integer
fAdj = (1.28 ^ 1.3 * N - 1) ^ 0.5
fRel = fL / f * fAdj
Bn = 0
Pn = 0
For k = 0 To N
c(k) = WorksheetFunction.Fact(2 * N - k) / _
(2 ^ (N - k) *
WorksheetFunction.Fact(k) * WorksheetFunction.Fact(N - k))
If WorksheetFunction.IsOdd(k) Then
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(Abs(k -
1) / 2, 0)) Then
Im(k) = -c(k)
Else
Im(k) = c(k)
End If
Pn = Pn + fRel ^ k * Im(k)
Else
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(k / 2, 0))
Then
Re(k) = -c(k)
Else
Re(k) = c(k)
End If
Bn = Bn + fRel ^ k * Re(k)
End If
Next k
BesselHPGain = 20 * WorksheetFunction.Log10(c(0)
/ (Bn ^ 2 + Pn ^ 2) ^ 0.5)
End Function
Public Function BesselBPGain(f As Double,
fL As Double, fU As Double, N As Double) As Double
Dim Im(0 To 25) As Double
Dim Re(0 To 25) As Double
Dim Bn As Double
Dim Pn As Double
Dim c(0 To 25) As Double
Dim fRel As Double
Dim fAdj As Double
Dim k As Integer
Dim fo As Double
Dim BW As Double
fAdj = (1.28 ^ 1.3 * N - 1) ^ 0.5
fo = (fU * fL) ^ 0.5
BW = (fU - fL) / fo
fRel = (((f / fo) - (fo / f)) / BW) * fAdj
Bn = 0
Pn = 0
For k = 0 To N
c(k) = WorksheetFunction.Fact(2 * N - k) / _
(2
^ (N - k) * WorksheetFunction.Fact(k) * WorksheetFunction.Fact(N - k))
If WorksheetFunction.IsOdd(k) Then
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(Abs(k
- 1) / 2, 0)) Then
Im(k) = -c(k)
Else
Im(k) = c(k)
End If
Pn = Pn + fRel ^ k * Im(k)
Else
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(k
/ 2, 0)) Then
Re(k) = -c(k)
Else
Re(k) = c(k)
End If
Bn = Bn + fRel ^ k * Re(k)
End If
Next k
BesselBPGain = 20 * WorksheetFunction.Log10(c(0)
/ (Bn ^ 2 + Pn ^ 2) ^ 0.5)
End Function
Public Function BesselBSGain(f As Double,
fL As Double, fU As Double, N As Double) As Double
Dim Im(0 To 25) As Double
Dim Re(0 To 25) As Double
Dim Bn As Double
Dim Pn As Double
Dim c(0 To 25) As Double
Dim fRel As Double
Dim fAdj As Double
Dim k As Integer
Dim fo As Double
Dim BW As Double
fAdj = (1.28 ^ 1.3 * N - 1) ^ 0.5
fo = (fU * fL) ^ 0.5
BW = (fU - fL) / fo
fRel = 1 / (((f / fo) - (fo / f)) / BW) * fAdj
Bn = 0
Pn = 0
For k = 0 To N
c(k) = WorksheetFunction.Fact(2 * N - k) / _
(2
^ (N - k) * WorksheetFunction.Fact(k) * WorksheetFunction.Fact(N - k))
If WorksheetFunction.IsOdd(k) Then
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(Abs(k
- 1) / 2, 0)) Then
Im(k) = -c(k)
Else
Im(k) = c(k)
End If
Pn = Pn + fRel ^ k * Im(k)
Else
If WorksheetFunction.IsOdd(WorksheetFunction.RoundDown(k / 2, 0))
Then
Re(k) = -c(k)
Else
Re(k) = c(k)
End If
Bn = Bn + fRel ^ k * Re(k)
End If
Next k
BesselBSGain = 20 * WorksheetFunction.Log10(c(0)
/ (Bn ^ 2 + Pn ^ 2) ^ 0.5)
End Function
Posted January 8, 2019
|
