Butterworth filter transfer function in time domain

Question

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:

Answer 1

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.

Answer 2

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}$$

Answer 3

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).

Answer 4

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;

}