5
$\begingroup$

My task is simple; I want to simulate analog low-pass filtering of an input signal, using Python. Note that the input signal is an array of values, not an analytical function.

My first question is if it simply is possible? If so, where is the line drawn between making a digital (FIR/IIR) filter and simulating an analog filter?

My second question is how I, in python, can make two equivalent simulations, using convolution in one, and converting the input signal to frequency domain in the other. Convolution is used in the link below, and it would help me a lot to see how exactly the same thing would be done by instead using fft (even though I think it's using an FIR filter, not simulating an analog filter).

http://glowingpython.blogspot.nl/2012/02/convolution-with-numpy.html

Thanks a lot in advance!

$\endgroup$
  • $\begingroup$ What do you know about the filter that is to be implemented? Impulse response, transfer function, other specifications? $\endgroup$ – Deve Mar 21 '13 at 8:29
  • 2
    $\begingroup$ Certainly I would draw the line at doing FFT based filtering. If you are simulating, you must want to preserve time domain characteristics and you won't do that with an FFT based approach. When you say “simulating” what are your simulation goals? Do you care about the filter structure? If so you'll need to model your simulation to preserve the filter structure as well as the time domain characteristics. So what is it that you really want to accomplish. That will determine what simulation method is appropriate. $\endgroup$ – user2718 Mar 22 '13 at 20:31
  • $\begingroup$ "Note that the input signal is an array of values" That sounds like digital filtering to me, not analog filtering. $\endgroup$ – endolith Mar 28 '13 at 19:50
  • $\begingroup$ Hi, Deve. I want to test different kinds of analog filters, but for example a 1st-order Bessel filter. $\endgroup$ – Uffe Mar 29 '13 at 13:56
  • 1
    $\begingroup$ @endolith. In my case it's not feasible to model the input with analytical functions, so I have no choice but to create it as an array. What's not clear to me is what the fundamental difference (if any) is between simulating an analog filter and making a digital filter. In scipy you can choose to make a filter analog or digital (see scipy.signal.bessel, so I guess there is a difference. $\endgroup$ – Uffe Mar 29 '13 at 14:35
7
$\begingroup$

What's not clear to me is what the fundamental difference (if any) is between simulating an analog filter and making a digital filter.

Either way, these functions will produce "ba" transfer function outputs, but the b and a are totally different.

For a 2nd-order filter, for instance, b = [b0, b1, b2] and a = [a0, a1, a2]. These are the coefficients of the transfer function.

For an analog filter, this represents a transfer function like this:

$$H(s) = \frac{b_2 s^2 + b_1 s + b_0}{a_2 s^2 + a_1 s + a_0}$$

For a digital filter, the transfer function is:

$$H(z)=\frac{b_0+b_1z^{-1}+b_2z^{-2}} {a_0+a_1z^{-1}+a_2z^{-2} }$$

Note that the input signal is an array of values

Then you are doing digital filtering. You can simulate analog filtering this way, but it's only going to be a simulation.

All the Python IIR filter functions you're talking about, when outputting a digital filter, design it as an approximation of the analog filter using an analog filter transfer function prototype and then converting it to digital using the bilinear transform and frequency warping. Yes, they can output analog transfer functions, but I think that's only useful for building actually analog electronic filters, or for plotting the frequency response, impulse response, etc.

The Bessel in particular is a bad approximation as you get near fs/4, because the important property of a Bessel filter is group delay, not amplitude, and phase/group delay is not preserved by the bilinear transform (red vertical line is fs/4):

enter image description here

More plots of digital Bessel filter

You can make the simulation closer to what an analog filter would produce by putting a large margin between the highest frequency in your signal and the sampling frequency, called oversampling. Assuming your "array of values" is correctly bandlimited, then you can upsample them and then apply the digital filter using lfilter and the ba digital filter coefficients you found earlier. If not, you need to bandlimit the analog signal and then sample it at a significantly higher rate. How did you get your sampled signal? If you are generating it digitally, it is probably not bandlimited, unless you specifically thought of bandlimiting while generating it.

for example a 1st-order Bessel filter

All 1st-order filters are the same, the "1st-order Bessel" is the same thing as a "1st-order Butterworth" or anything else.

$\endgroup$
  • $\begingroup$ Thanks for a very good answer! Currently I'm generating the input array within the python script, so I'll make sure to have significant oversampling. Also, currently I'm mostly interested in noise performance, so the phase response is of less importance. Currently I'm transforming the input to frequency domain, multiplying it with the filter frequency response and finally transforming back to time domain. Is that equivalent to using lfilter or filtfilt? Also, how could I use convolution to reach the same result? $\endgroup$ – Uffe Apr 2 '13 at 8:05
  • $\begingroup$ @Uffe: How are you generating the signal, though? Is it bandlimited? How are you calculating the frequency response? You can do filtering using frequency-domain multiplication, time-domain convolution, or using a direct filter implementation. All of these can be made to produce identical output, I believe, but there are implementation details. What are you trying to do? $\endgroup$ – endolith Apr 2 '13 at 14:38
  • $\begingroup$ Without going into too much detail, I'm generating something like a square wave but with changing amplitude levels. Then I add noise using 'random.normal'. I simply use fft to get the freq response. To convert back to time domain, however, I make sure to have hermitian symmetry. Finally I need to be able to reconstruct the original voltage levels. It seems to me that using fft and multiplying with the vector plotted to the left in your picture should give a very close estimation of applying an analog filter. This should avoid the fs/4 problem I guess. Do you agree? $\endgroup$ – Uffe Apr 4 '13 at 7:32
  • $\begingroup$ @Uffe: I'm guessing your "square wave" is not bandlimited? (switching between amplitudes is not good enough: music.mcgill.ca/~gary/307/week5/bandlimited.html and gist.github.com/endolith/407991) So when you generate it at a higher sampling rate, it will still not be bandlimited (square waves have infinite harmonics and your "oversampling" would just include more of them). Still, the higher in sampling frequency you go, the better approximation of analog you get, so I'd use lfilter and heavy oversampling. Increase the sampling rate until you don't see changes in output. $\endgroup$ – endolith Apr 4 '13 at 16:03
  • 2
    $\begingroup$ I think I understand what you mean. Any aliasing affects below the cutoff frequency of the filter I later apply will stay in the signal if I don't bandlimit my initial array. That is a difference compared to the analog situation which I'm trying to simulate. For simplicity I think I will use a high degree of oversampling. Thanks! $\endgroup$ – Uffe Apr 9 '13 at 12:03
1
$\begingroup$

Q1. Am I right in thinking that using scipy.signal.bessel with analog=1 and multiplying in the frequency domain will be the most similar to the analog situation?

What you get from scipy.signal.bessel with analog = 1 will be the continuous-time transfer function (assuming you want the ba version of the output). You will have to be careful to make sure you're doing the right sort of filtering: just using FFTs to simulate this will probably not work correctly.

I'm not an expert in python, but I suspect you want to use scipy.signal.freqs to get the continuous frequency response.

The trouble then is: how do you get the continuous frequency response of your input signal? Do you now how your analog input signal vector was generated?

EDIT

So reading the scipy manual, I see lsim or lsim2 are probably the options you want. It appears to allow you to apply a continuous-time system to data sampled at given sample points.

Q2. ... how [can] I, in python, can make two equivalent simulations, using convolution in one, and converting the input signal to frequency domain in the other

To do the convolution, you could just use scipy.signal.impulse to get the impulse response of your filter, and then use scipy.signal.convolve to convolve it with your sampled data.

To get the frequency domain version, you'll again need to understand how your analog input signal vector was generated.

$\endgroup$
  • $\begingroup$ Thanks for your answer. Regarding Q1, I can try to make my question specific. The picture in answer 1 above compares an analog freq response with a digital one. I assume that by using normal digital filtering (using e.g. lfilter) filtering will be made as seen in the right plots. If I instead transform to freq domain, and use the response seen to the left, I should get a closer approximation to an analog filter (assuming that I convert back to time domain correctly). Am I correct in this? $\endgroup$ – Uffe Apr 17 '13 at 11:55
  • $\begingroup$ Regarding Q2, also here to be specific, how can I make an exactly equivalent filtering, as seen in the link, but transforming to freq domain, and do the calculation there, instead of using convolution? (perhaps someone skilled in python can help here) $\endgroup$ – Uffe Apr 17 '13 at 12:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.