Simulating analog filter using convolution or converting with fft

Question

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!

Answer 1

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

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.

Answer 2

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.

Answer 3

I need something similar in Python but I can't understand what's wrong...

1. creating simple 1th order low-pass Butterworth filter, Fc=10Hz

# creating I order, 10 Hz, low pass Butterworth filter
fc = 10 # [Hz]
wc = 2*np.pi*fc # [rad/s]
b, a = signal.butter(1, wc, 'low', analog=True) # mi ritorna un sistema Num(s)/Den(s)
num=b
den=a

w, h = signal.freqs(b, a)   # ideal frequency response of filter
w = 0.5*w/np.pi             # [rad/s] to [Hz]
amp=20*np.log10(np.abs(h))  # amplitude to dB

'''
# ideal filter frequency response dB vs Hz
plt.plot(w, amp)
plt.xscale('log')
plt.grid()
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude [dB]")
plt.show()
'''

2. get impulse response

system=signal.lti(num, den)
t, impulse_resp = signal.impulse(system)
impulse_resp = impulse_resp[1:]     # removing first element due it's always 0.0
t = t[1:]

print("Risposta impulsiva:")        # filter impulse response (kernel of filter in time domain)
print(impulse_resp)

# normalizing kernel to 1.0
kernel=len(impulse_resp)
impulse_resp=impulse_resp /sum(impulse_resp)    # vado a normalizzare il kernel per il filtro
print(sum(impulse_resp))

3. creating white noise

# let's make some noise
mu, sigma = 0, 0.1
sample = 100000         # 100k samples
Fs=1000                 # [SPS] sampling frequency
x = np.arange(sample)
noise = np.random.normal(mu, sigma, len(x))  # noise

# time domain to frequency domain (noise)
#f, psd = signal.periodogram(noise, Fs)
f, psd = signal.welch(noise, Fs, scaling='spectrum')

# removing first element due related to 0 Hz (for log scale)
f=f[1:]
psd=psd[1:]

4. filtering in time domain using convolution and filter impulse response

# filtering in time domain using convolution and filter impulse response
filtrato=np.convolve(noise, impulse_resp)   # filtering noise
filtrato=filtrato[:len(filtrato)-len(impulse_resp)+1]   # equalizing length due to kernel dimension

# time domain to frequency domain for filtered noise
#f, psd2 = signal.periodogram(filtrato, Fs)
f, psd2 = signal.welch(filtrato, Fs, scaling='spectrum')

# as above
f=f[1:]
psd2=psd2[1:]

5. amplitude spectrum ratio to get filter frequency response

# ratio filtered/unfiltered signal
legenda=[]
tf=psd2/psd
plt.plot(f,20*np.log10(tf))
legenda.append("real")
plt.plot(w, amp)
legenda.append("ideal")
plt.xscale('log')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Attenuation [dB]')
plt.legend(legenda)
plt.grid()
plt.show()

Is this right or am I doing something wrong? Thanks