I am learning about the transfer function of IIR filters, and I am calculating it two different ways: directly from the formula and taking the z-transform of the output. They are giving different answers and I want to understand why. If X is the input and Y is the output, we have
$ |Y(z) |^2 = |H(z)|^2 |X(z)|^2$ so that
$ |Y(e^{i\omega T}) |^2 = |H(e^{i\omega T})|^2 |X(e^{i\omega T})|^2$. Does this imply that the power spectrum of the output is equal to the power spectrum of the input multiplied by $|H(e^{i\omega T})|^2$ ?
I tried to implement this with an actual signal and IIR filter. They differ but is it because of the initial ringing of the filter before it reaches a steady state? Is this initial ringing expected? Or am I missing a term somewhere? Otherwise is the method correct? When the input is a sine wave with a slightly different frequency, the green curve also looks like a gaussian but with a differing amplitude.
import numpy as np
from numpy.fft import rfft
from numpy import pi
from scipy.signal import lfilter
import matplotlib.pyplot as plt
b=[1,-1.414, 1.0]
a=[1,-1.273,0.81]
N_pts=1000
w=np.linspace(0,np.pi,N_pts//2)
z= np.exp(1j*w)
H = (b[0] + b[1]/z + b[2]/z**2)/(a[0] + a[1]/z + a[2]/z**2)
freq_vals=np.linspace(0,0.5,N_pts//2)
f=0.125 # frequency close to the notch
times = np.arange(N_pts)
signal = np.sin(2*np.pi*f*times)
filtered = lfilter(b,a,signal) # filtered signal
# note initial oscillations, then becomes zero as expected
## calculate |Y(z)|^2 two different ways
fft2_direct = np.absolute(rfft(filtered))**2
Yz_direct = fft2_direct[:-1] # ignore last freq
fft2_sig = np.absolute(rfft(signal))**2
fft2_sig = fft2_sig[:-1] # ignore last freq
Hmag = np.absolute(H)**2
Yz_H = fft2_sig*Hmag
plt.subplot(2,2,1)
plt.plot(freq_vals,np.absolute(H))
plt.title('|H(w)|^2')
plt.subplot(2,2,2)
plt.plot(filtered[:100])
plt.xlabel('first 100 time steps')
plt.title('filtered sine wave near notch')
plt.subplot(2,2,3)
plt.plot(freq_vals[1:],Yz_H[1:],'g-')
plt.plot(freq_vals[1:],Yz_direct[1:],'r--')
plt.title(' FFT')
plt.show()
