Supplement to my lectures on Linear Systems
R. Brigola, Georg-Simon-Ohm-Hochschule Nürnberg
Solution for the Exercise on the Design a Specified Butterworth Lowpass Filter
As an application example for a LTI-System design a Butterworth Filter (S. Butterworth, 1930):
At first, calculate the necessary order of the filter and the cutoff angular frequency according to its specification
as we have done in classroom.
( Reference for example: R. Brigola Fourieranalysis, Distributionen und Anwendungen, Vieweg, pg. 208 )
Specification:
To be concrete, calculate H for K=1, passband edge at 3 kHz, stopband edge at 5 kHz, minimal passband gain 0.9, maximal stopband gain 0.1.
Then, calculate the n zeroes , ... , with negative real parts of
, i.e .
Give the resulting frequency characteristics H of the Butterworth filter of order n and angular cutoff frequency (bandwidth).
Plot the frequency characteristics, the phase shift and delay of the filter.
H is given by
Solution:
At first calculate the necessary order of the filter and its cutoff frequency from the given specification
> | restart: |
> | K:=1.0; |
> | H1:=0.9; |
> | H2:=0.1; |
> | w1:=3000*(2*Pi); |
> | w2:=5000*(2*Pi); |
> | ord:=evalf(ln((K^2/H1^2-1)/(K^2/H2^2-1))/(2*ln(w1/w2))); |
> | n:=floor(ord)+1; # Order of the necessary Butterworth Filter |
> | Omega:=sqrt(w1*exp(-ln(K^2/H1^2-1)/(2*n))*w2*exp(-ln(K^2/H2^2-1)/(2*n))); |
> | evalf(Omega/(2*Pi)); # Cutoff frequency in Hz |
> |
H can be represented with the Butterworth polynomial in the form
The zeroes of Q with negative real parts are easily calculated by hand and give the zeroes of P
> | zero:=array(1..n); |
The zeroes of P are known to be ( you can compare them with the Maple solutions of Q(z)=0. Extract the zeroes with negative real part )
> | for k from 1 to n do |
> | zero[k]:=evalc(exp(I*((2*k-1)*Pi/(2*n)+Pi/2))); |
> | end do; |
Therefore the polynomial P is given by
> | P:=unapply(simplify(evalf(expand(product((s-zero[m]),m=1..n)))),s); |
> |
Now the frequency characteristics of the Butterworth lowpass is
> | H:=unapply(K/P(I*omega/Omega),omega); |
> |
> | plot(abs(H(t*Omega)),t=0..2,thickness=3,labels=["cutoff frequency scaled to 1",gain]); # plot the filter gain |
> | plot(abs(H(t*(2*Pi))),t=0..9000,colour=blue,thickness=3,labels=["frequency in Hz",gain]); # plot the filter gain |
The gain can also be represented as
> | H_abs:=unapply(K/sqrt(1+(w/evalf(Omega))^(2*n)),w); |
> | plot(H_abs(t*Omega),t=0..2,thickness=3,labels=["cutoff frequency scaled to 1",gain]); |
The frequency response can also be represented as a curve in the complex plane as follows.
Observe, that this is a qualitative graphics only unless you compute complex values for certain frequencies,
to see where you are on the curve for these.
> | with(plots):complexplot(H(w),w=0..20*Omega,thickness=3); |
We zoom the curve for high frequencies
> | complexplot(H(w),w=15*Omega..30*Omega,thickness=3); |
This shows approximately, that the curve near the origin crosses again the imaginary axis and approximates the origin from the left lower quadrant.
A mathematically well-educated engineer can see from that behavier of the curve - it begins for w=0 on the positive real axis,
goes in the mathematically negative direction to the origin, while crossing exactly as many quadrants of the complex plane as the filter order n (here n=6) indicates.
These facts tell that engineer, that the filter is stable (no resonance frequencies exist, a bounded input generates also a bounded output, it has the so-called BIBO-stability)
Remark for your future: A filter with rational frequency response is unstable, when the denominator polynomial has coefficients less or equal to zero !
Filters with frequency responses as above with respect to their start point for w=0, their orientation and number of quadrant crossings are stable.
Phase shift of the filter
> | assume(w,real); |
> | phi:=w->-sum(argument(I*w-Omega*zero[m]),m=1..n); |
> | plot(phi(t*Omega),t=0..2,colour=blue,thickness=3,labels=["cutoff frequency scaled to 1","Phase shift"]); |
> |
Group delay of the filter in seconds, you see that approximately the filter has only the superb delay of only 0.2 ms
> | delay:=w->-diff(phi(w*Omega),w); # observe, that the chain rule is used here, when the derivative is computed |
> | plot(delay(w)/Omega,w=0..2,colour=red,thickness=3,labels=["cutoff frequency scaled to 1","Group delay in s"]); |
> |
We observe that for frequencies < (1/2 * cutoff frequency) the gain and the delay are almost independent of the frequency.
Thus that filter provides high-fidelity transmission for an input with a bandwidth < Omega/2.
Let's test a few values
> | abs(H(Omega)); |
> | 20*log10(%); # damping in dB at the cutoff frequency |
> | abs(H(w1)); |
> | abs(H(w2)); |
> | 20*log10(%); # damping in dB at the stopband edge |
Let's calculate a table of the Butterworth polynomials P of order 2 to 10 as it can be found in references on circuit synthesis
> | restart: |
> | for n from 2 to 10 do |
> | for k from 1 to n do |
> | zero[k]:=evalc(exp(I*((2*k-1)*Pi/(2*n)+Pi/2))): |
> | end do; |
> | print('Order'=n); |
> | P:=simplify(evalf(expand(product((s-zero[m]),m=1..n)))): |
> | end do; |
> |
Learn more on filters and how filters can be realized by circuits. See for example
Wikipedia: Butterworth lowpass filters