f0=225kHz我想要一个带有和的数字 IIR 滤波器fs=53.125GHz。我可以想出传递函数并使用 Matlab 绘制它。当我尝试编写自己的代码来实现过滤器而不是 Matlab 的filter()函数时，就会出现问题。
第1步
%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)
到现在为止还挺好。接下来看一下时域响应
第2步
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')
正如预期的那样，有一条很长的尾巴。好的。现在让我们做一个 FFT 来交叉检查频率响应。
第 3 步
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')
看起来还是不错的。-3dB 点位于 224.57kHz。这告诉我，上面的时域响应是我想要得到想要的频域响应。
第4步
现在，当我尝试自己实现而不是使用filter(). 我采用传递函数并得出差分方程：
>> h_lpf
h_lpf =
1.331e-05 z + 1.331e-05
-----------------------
z - 1
Sample time: 1.8824e-11 seconds
Discrete-time transfer function.
看起来很简单，但这就是事情开始出错的地方
将分子和分母乘以 z^-1 得到因果形式：
1.331e-05 z + 1.331e-05 1.331e-05 + 1.331e-05 z^-1
----------------------- == -----------------------
z - 1 1 - z^-1
因此差分方程为：
y[n] = 1.331e-05*x[n] + 1.331e-05*x[n-1] + y[n-1]
直观地，我可以看到上面的方程实际上只是一个积分器。积分器不会产生所需的 LPF 响应。除了产生脉冲响应的衰减y[n-1]之外，还需要对另一个应用一些系数。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);
end
figure();
plot((0:200000-1)/fs, lpf_out(1:200000))
xlabel('Time [s]')
ylabel('Response')
title('Difference Equation Time Domain Response')
我知道我需要一个不同的系数y[n-1]。与分子系数相关的东西似乎很有希望。考虑 LPF（DC 增益为 0dB）将如何响应 DC 信号，我知道稳态输出将与输入相同：x[n]=y[n] 假设输入是一长串 1。2*1.331e-05由于 1.331e-05*x[n] + 1.331e-05*x[n-1]差分方程中的 ，每次过滤器迭代都在添加。y[n-1]除非缩放，否则输出将无限增加。
考虑输入只是 1 并且输出已经处于已知正确稳态的情况y[n-1] = 1：为了确保在y[n]=1存在1.331e-05*x[n] + 1.331e-05*x[n-1]前馈贡献的情况下输出保持不变，递归贡献必须是1-2*1.331e-05。因此，我们在这里结束：
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')
果然，这产生了与filter()步骤 3 中的方法相同的时域响应。但是，当我最初尝试从传递函数导出差分方程时，为什么没有到达这里？ 我是误解了传递函数还是误解了如何推导差分方程？