RealWorld Signals and Systems Case: Analog Filter Design with a Twist
You’re given the task of designing an analog (continuoustime) filter to meet the amplitude response specifications shown. You also need to find the filter step response, determine the value of the peak overshoot, and time where the peak overshoot occurs.
The objective of the filter design is for the frequency response magnitude in dB (20log_{10}H(f)) to pass through the unshaded region of the figure as frequency increases. The design requirements reduce to the passband and stopband critical frequencies f_{p} and f_{s} Hz and the passband and stopband attenuation levels A_{p} and A_{s} dB.
Additionally, the response characteristic is to be Butterworth, which means that the filter magnitude response and system function take this form:
Here, N is the filter order, f_{c} is the passband 3 dB cutoff frequency of the filter, and the poles, located on a semicircle is the lefthalf splane, are given by
This problem requires you to work in the frequency domain, the time domain, and perhaps the sdomain, depending on the solution approach you choose.
From the requirements, the filter frequency response has unity gain (0 dB) in the passband. The step response (a timedomain characterization) of the Butterworth filter is known to overshoot unity before finally settling to unity as
To design the filter, you can use one of two approaches:

Work a solution by hand, using the Butterworth magnitude frequency response H_{BU}(f) and the system function, H_{BU}(s).

Use the filter design capabilities of the SciPy signal package.
Finding the filter order and 3 dB cutoff frequency
Follow these steps to design the filter by using Python and SciPy to do the actual number crunching:

Find N and f_{c} to meet the magnitude response requirements.
Use the SciPy function N,wc=signal.buttord(wp,ws,Ap,As,analog= 1) and enter the filter design requirements, where wp and ws are the passband and stopband critical frequencies in rad/s and Ap and As are the passband and stopband attenuation levels (both sets of numbers come from the preceding figure). The function returns the filter order N and the cutoff frequency wc in rad/s.

Synthesize the filter — find the {b_{k}} and {a_{k}} coefficients of the LCC differential equation that realizes the desired system.
If finding circuit elements is the end game, you may go there immediately, using circuit synthesis formulas. Call the SciPy function b,a=signal.butter(N,wc,analog=1) with the filter order and the cutoff frequency, and it returns the filter coefficients in arrays b and a.

Find the step response in exact mathematical form or via simulation.
Here’s how to use the Python tools with the given design requirements and then check the work by plotting the frequency response as an overlay. Note: You can do the same thing in MATLAB with almost the same syntax.
In [379]: N,wc =signal.buttord(2*pi*1e3,2*pi*10e3,3.0, 50,analog=1) # find filter order N In [380]: N # filter order Out[380]: 3 In [381]: wc # cutoff freq in rad/s Out[381]: 9222.4701630595955 In [382]: b,a = signal.butter(N,wc,analog=1) # get coeffs. In [383]: b Out[383]: array([7.84407571e+11+0.j]) In [384]: real(a) Out[384]: array([1.00000000e+00, 1.84449403e+04, 1.70107912e+08,7.84407571e+11])
The results of Line [379] tell you that the required filter order is N = 3 and the required filter cutoff frequency is
The filter coefficient sets are also included in the results.
Use the real() function to safely display the real part of the coefficients array a because you know the coefficients are real. How? The poles, denominator roots of H_{BU}(s), are real or occur in complex conjugate pairs, ensuring that the denominator polynomial has real coefficients when multiplied out. The very small imaginary parts are due to numerical precision errors.
Checking the final design frequency response
To check the design, use the frequency response recipe.
In [386]: f = logspace(2,5,500) # log frequency axis In [387]: w,H = signal.freqs(b,a,2*pi*f) In [388]: semilogx(f,20*log10(abs(H)),’g’)
The figure shows the plot of the final design magnitude response along with the original design requirements.
Finding the step response from the filter coefficients
The most elegant approach to finding the step response from the filter coefficients is to find
The sdomain section of the figure tells you how to complete the partial fraction expansion (PFE) numerically. You have the coefficient arrays for H(s), so all you need to do is multiply the denominator polynomial by s. You can do this by hand or you can use a relationship between polynomial coefficients and sequence convolution.
When you multiply two polynomials, the coefficient arrays for each polynomial are convolved, as in sequence convolution.
Here, you work through the problem, using signal.convolve to perform polynomial multiplication in the denominator. To convince you that this really works, consider multiplication of the following two polynomials:
(x^{2} + x + 1)(x + 1) = x^{3} + 2x^{2} + 2x + 1
If you convolve the coefficients sets [1, 1, 1] and [1, 1] as arrays in Python, you get this output:
In [418]: signal.convolve([1,1,1],[1,1]) Out[418]: array([1, 2, 2, 1])
This agrees with the hand calculation. To find the PFE, plug the coefficients arrays b and convolve(a,[1,0]) into R,p,K = residue(b,a). The coefficients [1, 0] correspond to the sdomain polynomial s + 0.
In [420]: R,p,K = signal.residue(b,signal.convolve([1,0],a)) In [421]: R #(residues) scratch tiny numerical errors Out[421]: array([ 1.0000e+00 +2.3343e16j, # residue 0, imag part 0 1.0000e+00 +1.0695e15j, # residue 1, imag part 0 1.08935e15 5.7735e01j, # residue 2, real part 0 1.6081e15 +5.7735e01j]) # residue 3, real part 0 In [422]: p #(poles) Out[422]: array([ 0.0000 +0.0000e+00j, # pole 0 9222.4702 1.5454e12j, # pole 1, imag part 0 4611.2351 7.9869e+03j, # pole 2 4611.2351 +7.9869e+03j])# pole 3 In [423]: K #(from long division) Out[423]: array([ 0.+0.j]) # proper rational, so no terms
You have four poles: two real and one complex conjugate pair — a bit of a mess to work through, but it’s doable. Refer to the transform pair
to calculate the inverse transform for all four terms.
For the conjugate poles, the residues are also conjugates. This property always holds.
You can write the inverse transform of the conjugate pole terms as sines and cosines, using Euler’s formula and the cancellation of the imaginary parts in front of the cosine and real parts in front of the sine:
Putting it all together, you get y_{step}(t) = u(t) – e^{–9,222.47}^{t}u(t)–2 x 0.5774e^{–4,611.74}^{t}sin(7,986.89t)u(t). Having this form is nice, but you still need to find the function maximum for t > 0 and the maximum location. To do this, plot the function and observe the maximum.
A more direct approach is to use simulation via signal.lsim and the timedomain recipe. The system input is a step, so the simulation output will be the step response. From the simulated step response, you can calculate the peak overshoot numerically and see it in a plot. The IPython command line code is
In [425]: t = arange(0,0.002,1e6) # step less than smallest time constant In [426]: t,ys,x_state = signal.lsim((b,a),ones(len(t)),t) In [428]: plot(t*1e3,ys)
Using the time array t and the step response array ys, you can use the max() and find() functions to complete the task:
In [436]: max(real(ys)) # real to clear num. errors Out[436]: 1.0814651457627822 # peak overshoot is8.14% In [437]: find(real(ys)== max(real(ys))) Out[437]: array([534]) # find peak to be at index 534 In [439]: t[534]*1e3 # time at index 534 in ms Out[439]: 0.5339