2

I have a 1st order low-pass filter (LPF) in the frequency domain and I want to digitize it. I am comparing the frequency response graphs for testing, but I'm getting weird results...

Though very basic, I couldn't get it right from reading the scipy.signal.bilinear help page or the net.

  • num, den: numerator and denominator in S plane
  • b, a: I'm expecting it to be b, a coefficients of a digital difference equation filter (IIR), of the shape: Y[n] = a0*X[n] + a1*X[N-1] + ... - b1*Y[n-1] ...

Here is a code example:

Fs = 48000.0
f = 2 * np.logspace(1,4,1024)

num =  [0 , 1]
den = [0.001 , 1]

tmp, H = sig.freqs(num, den, worN=1024)
b, a = sig.bilinear(num, den, 1.0)
tmp, Hd = sig.freqz(b,a, worN=1024)

plt.semilogx(f, 20*np.log10(np.abs(H)))
plt.semilogx(f, 20*np.log10(np.abs(Hd)))

What am I doing wrong?

Jongware
  • 22,200
  • 8
  • 54
  • 100
Malcolm Rest
  • 103
  • 7

2 Answers2

1

The problem is that you are not using the tmp in the x axis when plotting, and also the freqz gives the normalized tmp vector in radians/sample:

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

Fs = 48000
num =  [0 , 1000]
den = [1 , 1000]

w1, H = sig.freqs(num, den, worN=1024)
b, a = sig.bilinear(num, den, Fs)
w2, Hd = sig.freqz(b, a, worN=1024)

fig = plt.figure()
plt.title('Filter frequency response')
plt.semilogx(w1, 20*np.log10(np.abs(H)),'b')
plt.semilogx(w2*Fs, 20*np.log10(np.abs(Hd)),'k')
plt.ylabel('magnitude [dB]')
plt.xlabel('frequency [Hz]')
plt.grid()
plt.axis('tight')
plt.xlim([0.001, Fs/2])
plt.show()

This code works perfectly fine. Hope it helps.

SalvadorViramontes
  • 540
  • 4
  • 18
0

Yes, this is very helpful. I then took it to the next step so that both 'w' vectors are the same lenth and cover all of the interesting sampling points, like so:

Fs = 48000.0
f = 2 * np.logspace(1,4,1024)
w = 2 * np.pi * f

Snum =  [0 , 1]
Sden = [0.001 , 1]

w1, H = sig.freqs(Snum, Sden, worN=w)
b, a = sig.bilinear(Snum, Sden, Fs)
w2, Hd = sig.freqz(b,a, worN=w1/Fs)

Also, if the prsented graph is in Hz (not Rad/Sec0), w should be divied by 2*PI:

plt.semilogx(w1/np.pi/2, 20*np.log10(np.abs(H)), 'r')
Malcolm Rest
  • 103
  • 7