# Preliminaries

I need to know the response of Butterworth filter in time domain of the following input (at least numerically): $$x(t) = \theta(t)(1-\exp\{-at\})\exp\{-bt\}, \quad a,b>0,\quad a\gg b$$ In order to find the response I decided to find the transfer of Butterworth filter in time domain and then convolve it with $$x(t)$$ numerically. But the result I've got doesn't seem to be correct. Here how it went.

# Transfer function in $$s$$-domain

Transfer function of Butterworth filter is the following: $$H(s) = \prod_{k=1}^{n}\frac{1}{(s-s_{k})}, \quad s_{k} = \omega_{c}e^{i\frac{2k+n-1}{2n}\pi}$$ So in order to find transfer function in time domain I used this equality: $$h(t) = \sum_{s_{k}}\mathrm{res}H(s_k)\exp\{s_k t\}$$

# Transfer function in $$t$$-domain

First of all, I found that $$\Im{h(t)} = 0$$ i.e. it is real. I used this fact as a cross check. I wrote the code that calculates $$h(t)$$ using complex numbers. I found the following result: The question is : could somebody verify if it looks right? I haven't found any pictures of the Butterworth transfer function in time domain.

# P.S

I could post the code I used to produce those results, but keep in mind it is based on ROOT CERN. You can see the code and tests here.

# Update

I have performed one more test: calculation of the magnitude of the frequency response of Butterworth filter of the order $$n=3$$ and cut frequency $$\omega_{c}=2$$ rad/s. The result is perfect: • I updated my answer to include the h(t) expressions for M = 2 and 3. Maybe this will help you.
– Ed V
Feb 11, 2020 at 1:09

A plot of the normalized impulse responses, for the n = 2 through 10 Butterworth low pass filters, are given by H.J. Blinchikoff, A.I. Zverev, "Filtering in the Time and Frequency Domains", Wiley-Interscience, John Wiley & Sons, NY, ©1976, p. 113. This is shown below. They do not give the h(t) expressions in the book, at least where I have looked thus far. They say the figure is from Zverev's earlier book: A.I. Zverev, "Handbook of Filter Synthesis", Wiley, NY, ©1967. EDIT: Looking back in my old computation notebooks, I found the impulse expression for the 3rd order Butterworth low pass filter. This pic is from 1983: Hopefully, this is algebraically the same as your result.

EDIT 2: Gregg 1 gives $$h(t)$$ for the $$M = 2$$ Butterworth LPF: $$h(t) = e^{-x}sin(x)U(t)$$ where $$x = t/τ√2$$, $$τ = RC = 1/(2πB)$$ and $$B$$ is the corner frequency in Hz. Butterworth LPFs are fine for many purposes, but not particularly good for others, e.g., they have relatively poor 1% settling times with step inputs 2.

References:

1. W.D. Gregg, “Analog and Digital Communication”, Wiley, NY, ©1977.

2. E. Voigtman, J.D. Winefordner, “Low-pass filters for signal averaging”, Rev. Sci. Instrum. 57 (1986) 957-966.

EDIT 3: Even though Matt has provided an excellent answer, for completeness I decided to test my 1983 $$h(t)$$ expression (see above) by comparing it with the numerically generated impulse response from a computer simulation $$H(s)$$ transfer function block. As per the OP, I assumed $$2πB = 1$$. The next figure shows that the results are the same. • Thank you for your answer. I'll find the literature you refer here, but at the first sight the graph here is identical to one I've made. I' will put more details and additional results in my question then. Feb 10, 2020 at 15:45
• Ok. I have checked the case $M=2$. And it is exactly what I get numerically with my algorithm except it should be normalized by $\sqrt 2$ multiplier, i.e. $h(t) = \sqrt 2 e^{-x}\sin{x}$. Feb 11, 2020 at 4:24
• Of course, I assume $2\pi B = 1$. Feb 11, 2020 at 5:08
• I just did the $M = 2$ case, like the $M = 3$ case above, with $2πB = 1$, i.e., $τ = 1$. The results confirm that the normalizing factor is $√2/τ$, which is simply $√2$ when $2πB = 1 = τ$. So either Gregg forgot the normalizing factor (it was a library book I do not have) or I did. This has been fun, I corrected a 37 year old error and thanks for posting the question!
– Ed V
Feb 11, 2020 at 23:08
• I also have compared your formula for $M=3$ with my algorithm. And it is fits perfectly. See P.S. section in the post. Feb 12, 2020 at 7:37

There's no need to use numerical methods here. The most straightforward way to compute the output is to see that the filter's impulse response is given by

$$h(t)=\sum_{k=1}^Nr_ke^{s_kt}u(t)=\sum_{k=1}^Nh_k(t)\tag{1}$$

where $$N$$ is the filter order, $$u(t)$$ is the unit step function, and $$r_k$$ are the coefficients of the partial fraction expansion of $$H(s)$$:

$$H(s)=\sum_{k=1}^N\frac{r_k}{s-s_k}\tag{2}$$

If you have an input signal of the form $$x(t)=e^{\alpha t}u(t)$$ - and your input signal is just the sum of two such signals with different exponents - then the output signal can be written as

$$y(t)=(x\star h)(t)=\sum_{k=1}^N(x\star h_k)(t)\tag{3}$$

And the convolutions $$(x\star h_k)(t)$$ can be easily computed analytically:

\begin{align}(x\star h_k)(t)&=r_k\int_{-\infty}^{\infty}e^{\alpha \tau}u(\tau)e^{s_k(t-\tau)}u(t-\tau)d\tau\\&=r_ke^{s_kt}u(t)\int_{0}^{t}e^{(\alpha-s_k)\tau}d\tau\\&=r_k\frac{e^{\alpha t}-e^{s_kt}}{\alpha-s_k}u(t),\qquad \alpha\neq s_k\tag{4}\end{align}

• This is what I was looking for. I guess OP, too. Feb 11, 2020 at 10:32
– Ed V
Feb 11, 2020 at 13:25
• Yes, of course, it is possible to find $h(t)$ analytically, but I think it's a pain to find all $r_k$'s unless you know the painless way to calculate them for arbitrary $n$. That is why I chose to find it numerically. Feb 12, 2020 at 4:35
• @LRDPRDX: The analytic solution is really simple in this case, so I don't see why one would want to put up with inaccuracies due to numerical approximations. The residues are computed by evaluating polynomials. Have a look at Matlab's/Octave's residue function. In your case it's even easier because you know that there are only single poles. Feb 12, 2020 at 10:40

If you want to have a "numerical grasp" and you're not afraid of getting a little bit dirty, you can check the numbers with LTspice. I don't know how well you know to work with it, so I'll just explain it, feel free to ignore all the redundant info.

Here you can download the archive, out of which you only need Filter.asy and filter.sub. Create a new schematic and save it in a folder of your choice. Then copy the two files in that same directory as the location of the saved schematic. Then, with the schematic opened, press F2 and, from the Top Directory drop-down menu, select the current folder's location. The Filter symbol should be available to be placed now. From the same drop-down menu, return to the default location, where you can select voltage, place it next to the filter. Press g to place the ground symbols and F3 to add the wires, as seen in the picture: Right-click on the Filter and double-click on Bessel to bring up a drop-down menu, from which you can select Butterworth (not Butterorth_IIR). Set fp1=1, fs2 and fp2 both to zero, and fs1>fp1 (any value). On the last line, set N={x}. The filter is set now, rename the source as seen in the picture (right-click on the V), press s and add the .step text (left-click to place), then again with the .tran text. You should be set to run (the running man icon in the toolbar). After the simulation, with left-click on the desired node you can plot the signal.

Then you can left-click on the label of the trace in the waveform window and that will bring up a cursor that can be moved and you can read all the numbers you need. To place the cursor on different traces, use the up/down arrows. If you only need one trace, ctrl-right-click on the .step command and select comment, then right-click on the Filter and set N=<...>. For more usage you can read some here (not sure if it's readable/understandable enough).

If you want to go this way, good luck. :-)

PS: Forgot to say that the corner frequency is set by fp1, so if you want in radians, just set it to fp1={1/(2*pi)} (and fs1>fp1).

• Could I also find the response of my input (see post)? Feb 11, 2020 at 4:29
• Sure. Instead of a voltage source, add a behavioural source (F2 > bv) and add this expression: (1-exp(-a*time))*exp(-b*time), then press s and add .param a=<...> b=<...>. These can be changed at any time. If that $\theta(t)$ is the step function, you're done, if it's a sine, add sin(2*pi*f*time) (replace f with your frequency), or any other function. For extra precision, press s and add .opt plotwinsize=0, and impose a small enough timestep by right-clicking on the .tran command (a popup will appear). You can go as high as N=32 and it can be all/low/high/bandpass, bandstop. Feb 11, 2020 at 8:13
• I know that I have seen it somewhere, derived, the step response for generic transfer function, and it involved decomposing the t.f. into partial fractions as a sum, it was all parametric, based on the roots of the t.f.. I can't find it now, but I'll search some more. Could you post a link to your code, pastebin/etc? Feb 11, 2020 at 8:35
• @LRDPRDX Matt's answer is what you're looking for, but I'd still like to see the code, if you think you can post it. Feb 11, 2020 at 10:44
• @a concerned citizen, see the P.S. section in the post. I added a link to github there. Feb 12, 2020 at 7:38

Hi, a bit late to the game, but a solution for an Nth-order Low-pass Butterworth time-domain response written in C# is shown below. I hope it will (still) be of use.

double ButterLowPass(double time, double Fmax, int order)


{ double angle, cosPart, expPart, h_of_t, t, scale;

Complex Delta, RootN, Roots, Residue;  // Best to use complex numbers, to work with the roots
int ir, ii;

h_of_t = 0.0;
scale = 2.0 * Math.PI * Fmax;

// If time isn't normalized, the filter corner point would sit at w =  1 rad/s
t = time * scale;

for (ir = 1; ir <= order / 2; ir++)
{
angle = Math.PI * (order + (2.0 * ir) - 1.0) / (2.0 * order);
RootN = Complex.FromPolarCoordinates(1.0, angle);
Roots = Complex.One;

for (ii = 1; ii <= order; ii++)
{
if (ii == ir) continue;

angle = Math.PI * (order + (2.0 * ii) - 1.0) / (2.0 * order);

Delta  = RootN - Complex.FromPolarCoordinates(1.0, angle);
Roots *= Delta;

}

// after multiplying all poles, we need to put them in the denominator
Residue = Complex.One / Roots;
cosPart = 2.0 * Residue.Magnitude * Math.Exp(RootN.Real * t) * Math.Cos(RootN.Imaginary * t + Residue.Phase);
h_of_t += cosPart;
}

if ((order % 2) == 1)
{
// compute the real exponential component for odd - order responses

Roots = Complex.One;

ir = (order + 1) / 2;
angle = Math.PI * (order + (2.0 * ir) - 1.0) / (2.0 * order);
RootN = Complex.FromPolarCoordinates(1.0, angle);

for (ii = 1; ii <= order; ii++)
{
if (ii == ir) continue;

angle = Math.PI * (order + (2.0 * ii) - 1.0) / (2.0 * order);

Delta  = RootN - Complex.FromPolarCoordinates(1.0, angle);
Roots *= Delta;

}

// after multiplying all poles, we need to put them in the denominator
Residue = Complex.One / Roots;
expPart = Residue.Real * Math.Exp(-t);
h_of_t += expPart;

}

return h_of_t * scale;


}

• Sorry, it seems that formatting has been messed up a bit. Not sure why... Feb 24, 2021 at 21:14