How to Use Filter Equations in Software and Spreadsheets |
||
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 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 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 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 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 Related Pages on RF Cafe - Butterworth Lowpass Filter Gain, Phase, and Group Delay Equations - Butterworth Highpass, Bandpass, & Bandstop Filter Gain, Phase, and Group Delay Equations - How to Use Filter Equations in a Spreadsheet - Filter Equivalent Noise Bandwidth - Filter Prototype Denormalization - Bessel Filter Prototype Element Values - Butterworth Lowpass Filter Poles - Butterworth Filter Prototype Element Values - Chebyshev Lowpass Filter Poles - Chebyshev Filter Prototype Element Values - Monolithic Ceramic Block Combline Bandpass Filters Design - Coupled Microstrip Filters: Simple Methodologies for Improved Characteristics |
||