I want a digital IIR filter with
fs=53.125GHz. I can come up with the transfer function and plot it using Matlab. The problem arises when I try to do my own code that implements the filter instead of Matlab's
%Define transfer function lpf_f0 = 225e3; fs = 53.125e9; tau = 1/(2*pi*lpf_f0 ); lpf_num = 1; lpf_den = [tau 1]; %CT-DT conversion [lpf_numd,lpf_dend] = bilinear(lpf_num, lpf_den, fs); %Bode Plot h_lpf= tf(lpf_numd,lpf_dend, 1/fs); bopts = bodeoptions; bopts.FreqUnits = 'Hz'; bopts.PhaseVisible = 'off'; bode(h_lpf, bopts)
So far so good. Next, take a look at the time-domain response
data = zeros(2^30,1); data(100) = 1; lpf_out_filt = filter(lpf_numd,lpf_dend, data); figure(); plot((0:200000-1)/fs, lpf_out_filt(1:200000)) xlabel('Time [s]') ylabel('Response') title('LPF Time Domain Response')
As expected, there's a really long tail. Good. Now let's do an FFT to cross-check the frequency response.
lpf_out_cplx = fftshift(fft(lpf_out_filt)); f_step = fs/length(lpf_out_cplx); lpf_out_cplx = lpf_out_cplx(length(lpf_out_cplx)/2:end); lpf_out_mag = abs(lpf_out_cplx); lpf_out_freq = (1:length(lpf_out_mag))*f_step; lpf_out_mag = 20*log10(lpf_out_mag); i_f0 = find(lpf_out_mag<=-3,1); figure(); semilogx(lpf_out_freq(1:i_f0*10), lpf_out_mag(1:i_f0*10)) title('LPF Response') xlabel('Frequency [Hz]') ylabel('Magnitide (dB)') text(lpf_out_freq(i_f0), lpf_out_mag(i_f0), 'x')
Still looking good. -3dB point is at 224.57kHz. This tells me that the above time-domain response is what I want to get the desired frequency-domain response.
Now, the trouble comes when I try to do my own implementation instead of using
filter(). I take the transfer function and come up with the difference equation:
>> h_lpf h_lpf = 1.331e-05 z + 1.331e-05 ----------------------- z - 1 Sample time: 1.8824e-11 seconds Discrete-time transfer function.
Seems straighforward, but this is where things start to to awry
Multiply by numerator and denominator by z^-1 to get it in a causal form:
1.331e-05 z + 1.331e-05 1.331e-05 + 1.331e-05 z^-1 ----------------------- == ----------------------- z - 1 1 - z^-1
Difference equation is thus:
y[n] = 1.331e-05*x[n] + 1.331e-05*x[n-1] + y[n-1]
Intuitively, I can see that above equation is really just an integrator. An integrator won't yield the desired LPF response. There needs to be some coefficient applied to the
y[n-1] other than
1 to yield the decay in the impulse response. Sure enough, I get a step-like response to a input pulse:
lpf_out = zeros(size(data)); lpf_out(1) = data(1)*lpf_numd(1); for i=2:length(data) lpf_out(i) = lpf_numd*data(i-1:i) + lpf_out(i-1); end figure(); plot((0:200000-1)/fs, lpf_out(1:200000)) xlabel('Time [s]') ylabel('Response') title('Difference Equation Time Domain Response')
I know I need a different coefficient on the
y[n-1]. Something related to the numerator coefficients seems pomising. Thinking about how a LPF (with 0dB gain at DC) will respond to a DC signal, I know the stead-state output will be the same as the input:
x[n]=y[n] Assuming the input is a long string of 1's. Every filter iteration is adding
2*1.331e-05 on account of the
1.331e-05*x[n] + 1.331e-05*x[n-1] in the difference equation. The output will increase unbounded unless the
y[n-1] is scaled.
Consider the case where the input is just 1's and the output is already at the known correct stead-state:
y[n-1] = 1. To make sure the output stays at
y[n]=1 in the presence of the
1.331e-05*x[n] + 1.331e-05*x[n-1] feed-forward contriubtion, the recursive contribution must be
1-2*1.331e-05. Thus we end up here:
y[n] = 1.331e-05*x[n] + 1.331e-05*x[n-1] + (1-2*1.331e-05)*y[n-1]
lpf_out = zeros(size(data)); lpf_out(1) = data(1)*lpf_numd(1); for i=2:length(data) lpf_out(i) = lpf_numd*data(i-1:i) + lpf_out(i-1)*(1-sum(lpf_numd)); end figure(); plot((0:200000-1)/fs, lpf_out(1:200000)) xlabel('Time [s]') ylabel('Response') title('Difference Equation Time Domain Response')
Sure enough, this yeilds the same time-domain response that the
filter() approach in Step 3. However, why didn't I get here when I initially tried to derive the difference equation from the transfer function? Was I misinterpreting the transfer function or misunderstanding how to derive the difference equation?