【模拟滤波与数字滤波区别】--双线性变换的预畸变处理理解及验证；matlab中freqz、freqs、loglog、bode函数绘图的区别

前言

作者耗时4天，翻阅了不低于20篇文章，关于在模拟滤波及数字滤波设计时遇到的问题，一点点总结出相应的解决及验证方法；希望能帮助广大读者，节约时间。
主要问题包括如下：

• 双线性变换为什么要进行预畸变处理，预畸变处理后为什么我绘制的bode图截止频率对不上？[^1]
• 为什么系统使用c2d函数或者d2c函数后，绘制的bode图不一致？
• 为什么loglog和plot函数的图特么对不上，不理解？
  最后给出了matlab测试代码，包含如下内容：
• freqz、freqs、loglog、bode、plot函数在数字滤波和模拟滤波中的使用及区别（有注释）
• 基于低通滤波器设计 双线性变化进行预畸变与不进行预畸变处理和区别。
• 基于带通滤波器设计 双线性变化进行预畸变与不进行预畸变处理和区别。

解决方法及验证代码展示

第一个问题

• 首先说模拟滤波器和数字滤波器的区别
  模拟滤波器针对连续系统，一阶低通滤波器s域表达式$\frac{10}{s+10}$，就是模拟滤波器；将其进行z变换的表达式$\frac{0.004975z+0.004975}{z-0.99}$，就是对应的数字滤波器，针对离散系统。我们实际的采样数据都是离散的，一般用都是数字滤波，再者这里z变换表达式式是受采样频率影响的，z变化有2种方式，这是接下来要说的。
• z变换的方式
  上边说了s域转z域需要z变换，z变换一般有如下多种方式：
  
  对于冲激响应不变法和双线性不变化法这两种方法参考原理参考。我这里只介绍双线性变换法。以下是双线性变换法的一个优缺点
  
• 双线性变换为什么要进行预畸变处理
  双线性变换法存在严重的非线性频率变换，会发生畸变。怎么理解这个畸变，举个栗子，期望滤波器的截止频率$w_d$，这是用于实际采样数据的，所以是一个数字滤波器的截止频率，而滤波器设计都是从s域开始再转z域的这么一个过程，如果直接用$w_d$去设计s域的模拟滤波器，模拟滤波的截止频率是$w_d$，没问题，但是进行双性变换转z域的数字滤波，截止频率就发生了变化，不在$w_d$处了，这就是畸变。
• 如何进行预畸变处理？
  简单，但是这里要注意频率(Hz)与角频率(rad/s)的区别。我们的目的是使设计出来的数字滤波器截止频率为$w_d$。那么预畸变的模拟滤波器的截止频率（角频率rad/s）$w_a=\frac{2}{T}tan\frac{w_{d}T}{2}$。而对应的基于频率（Hz）的转换关系就是除以$2\pi$，即$w_{ahz}=\frac{2}{T\ast 2\pi }tan\frac{w_{dhz}\ast 2\pi \ast T}{2}=\frac{1}{T\pi }tan(w_{dhz}\pi T)$。
  然后有些人和我一样对上图有疑惑，为什么别人写的是$=\frac{2}{T}tan(\frac{w}{2})$，那是因为他这里w指的就是模拟滤波的截止频率经双线性变化压缩后的模拟频率，这个是$=w_{d}T$的。所以你放心，结果是没有错的，如果你对原理很有执念，参考。
• 预畸变处理后为什么我绘制的bode图截止频率对不上？
  那么你肯定是使用c2d函数或者绘图的时候出现了问题，这就涉及到了第二各大问题

第二个问题

• 为什么系统使用c2d函数或者d2c函数后，绘制的bode图不一致？
  因为在使用基于’tustin’方法转换时，发生了畸变，所以会不一样，且基于双线性变换的离散系统bode图的频率范围为$\pm \pi /T$,绘制出来会有一根竖线，如
  
  然后一个重要的点是，对于没系统学过的朋友，可能对freqs/freqz/bode/loglog这些matlab函数比较陌生，网上也没人比较全面的介绍，这可能让你分析幅频图像时遇到一些困难，我也是一样。那我告诉你，这些函数绘制出来的图主要区别在于横纵轴坐标的刻度、单位有所区别。这些可以在我的测试代码中看到。

第三个问题

同第二问题一样，这些函数的绘图区别在于横纵轴坐标的刻度、单位，这些可以在我的测试代码中看到。

matlab测试代码展示

给出测试的matlab代码，涉及的问题及答案都在代码中。

%% 滤波器练习
clear;clc;clf;
% wp = 2;ws = 4;
% rp = 1;rs = 30;
% fs = 100;
% 
% %传入的参数分别为：通带、阻带截止频率；通带、阻带衰减系数
% %返回参数为：阶数、3dB截止频率
% [N,wc] = buttord(wp,ws,rp,rs,"s");
% Wc = [wc,wc+1.5];
% 
% %计算多项式形式的分母分子系数
% [B1,A1] = butter(2,wc,'s'); %模拟低通滤波器
% [B2,A2] = butter(2,wc,"high","s");%模拟高通滤波器
% [B3,A3] = butter(1,Wc,'stop','s');%模拟带阻滤波器 对于带通和带阻设计，n 表示滤波器阶数的一半。
% [B4,A4] = butter(N,Wc,'s');%模拟带通滤波器（最后一个参数's'表示模拟，去掉后为数字为z域）
% 
% w = 0:0.001:6;

% for i=1:4
%     %创建变量名
%     var_name = strcat('H',num2str(i));
% 
%     %计算结果
%     res = freqs(eval(['B',num2str(i)]),eval(['A',num2str(i)]),w);
%     %取对数数
%     res = 20*log10(abs(res));
%     %将计算结果放到变量名中
%     eval([var_name,'=res;']);
% 
%     figure(i);
%     plot(w,eval(['H',num2str(i)]));
%     grid on;
%     xlabel('频率/(rad/s)');ylabel('幅度/dB');pow2(3,3)
% end


% [b,a]=bilinear(B3,A3,1000);
% 
% tf_fun = tf(b,a);
% figure(5)
% bode(tf_fun,'b');hold on;grid on;

%% 数字滤波和模拟滤波绘制的幅频曲线图区别及分析
% *********************一、从数字滤波的角度去分析，freqz与bode绘制的图像区别**************************************
% 结论：1.对于离散系统，根据双线性变换原理，其频率范围为 ±π/T。
%      2.butter在设计离散滤波器的时候进行了预畸处理，所以其离散系统bode图截止频率在截止频率30Hz处，同时对应的连续系统Wc不在30Hz处
%        角频率转换关系：Ω = 2/T *tan(W*T/2)。所以当W*T/2足够小时，Ω=W，为线性转换关系，基于bilinear转换的离散和连续系统的失真比
%        较小；另一个角度可以看出T过大或者W过大，都会加剧失真
%      3.freqz得到的频率W1需要进行转换，绘制的图才能和bode绘制图对应上(数学关系暂未推导) 
%      4.绘图注意：刻度和单位区分开；图1，2，3的区别：直接freqz绘制的图1与图3一致；图2 xy刻度为指数10^k；图3 x刻度为k,y刻度为db;
%        图4：标准bode图(x刻度10^k，y刻度db)
%      5.从图4可以看出离散系统用tustin转连续系统的bode图其截止频率是畸变了的。离散的为30*2pi，连续系统为275 rad/s。
fs = 50;
w_1 = 0:0.1:100;
% 可以看到butter函数源码是使用了双线性变化，并对截止频率进行预畸处理
[b1,a1] = butter(2,30/(fs/2)); %数字滤波DF，截止频率为30Hz; 
sys_1 = tf(b1,a1,1/fs);
[H1,W1]=freqz(b1,a1);
mag1=abs(H1); % 取幅度值实部
pha1=angle(H1); % 取相位角
db1=20*log10(mag1+eps); % 转换为分贝
f1=W1*fs/(2*pi); % 将 数字角频率 转为对应的 模拟频率Hz（w=ΩT*2pi-> Ω=w*fs/(2pi)）

figure(1)
freqz(b1,a1,w_1,fs);hold on;
figure(2)
loglog(f1,mag1);
figure(3)
plot(f1,db1);
figure(4)
sys_1_2c = d2c(sys_1,'tustin');
% bode图绘制的是以角频率为横坐标的单位，为Hz*2*pi
bode(sys_1,'r',sys_1_2c,'b');hold on;legend('离散系统','连续系统');
% ********************************************************************************************************


% *****************二、从模拟滤波的角度去分析，freqs与bode绘制的图像区别**************************************
% 结论: 1.绘图注意：freqz和freqs绘制出的图刻度不一致需区别开，freqz刻度是：x为k,y为db；freqs刻度是：xy都是10^k
%         其他同上述离散系统的分析一样；
%       2. loglog(w_1,mag2) = plot(w_1,20*log10(mag2)),区别是横坐标的刻度，前者是10^k一个刻度，后者k一个刻度
[b2,a2] = butter(2,40,'s'); %模拟滤波AF,截止频率为40rad/s
sys_2 = tf(b2,a2);
w_2 = logspace(-1, 3);
[H2,W2]=freqs(b2,a2,w_2); % W2为角频率rad/s
mag2=abs(H2);% 取幅度值实部
pha2=angle(H2);% 取相位角

figure(1)
freqs(b2,a2,w_2)% 绘制的是幅值对数图，分子分母刻度都是10^
figure(2)
loglog(w_2,mag2);% 分子分母刻度都是10^k
figure(3)
plot(w_2,20*log10(mag2))% 分子分母刻度都是均匀刻度k，所以与图1有区别，但是各点的值都对得上
figure(4)
bode(sys_2);hold on; % 绘制的是标准bode图(x刻度10^k，y刻度db)
% ***********************************************************************************************************

%% 双线性预畸变处理前后区别及分析
% ***********************************一、双线性预畸变处理频域验证基于低通****************************************
%！ ** matlab对sys_2进行双线性变换变成离散系统 **
% 结论：1.图5中离散系统相对连续系统依然发生了截止频率畸变，和前面分析一样，且离散的频率<连续的频率。离散映射到连续的频率
%       为 wa=2/T*tan(wd*T/s)。
%       2.sys_2_c2d和sys_3都是sys_2进行 s=2/T * (1-z^-1)/(1+z^-1)得来的，两者相等，验证了c2d函数的实现。不涉及预畸变
%       3.图6可以得到预畸变处理的离散系统bode图符合预期的40rad/s的截止频率
Ts = 1/fs;
sys_2_c2d = c2d(sys_2, Ts, 'tustin');
figure(5)
bode(sys_2_c2d,'r--',sys_2,'b');hold on; legend('离散系统','连续系统');
%！ ** 手动对sys_2进行双线性变换变成离散系统,不对wc预畸变处理 **
yinzi =(4/Ts^2) + 2*56.57/Ts+1600;
aa1 = 1;
aa2 = (-8/Ts^2 + 3200)/yinzi;
aa3 = ((4/Ts^2)-2*56.57/Ts+1600)/yinzi;
bb1 = 1600/yinzi;
bb2 = 3200/yinzi;
bb3 = 1600/yinzi;
sys_3 = tf([bb1 bb2 bb3], [aa1 aa2 aa3], Ts);

%！ **对sys_2进行双线性变换变成离散系统,预畸变处理,希望数字滤波的截止频率为40 rad/s**
wc_1 = 2/Ts * tan(40*Ts/2);
[b6, a6] = butter(2,wc_1, 's');
sys_6_s = tf(b6,a6);
bode();hold on;
sys_6_z = c2d(sys_6_s, Ts, 'tustin');
figure(6)
bode(sys_6_z,'r--',sys_6_s,'b');hold on; legend('离散系统','连续系统');
% ***********************************************************************************************************


% *****************************二、双线性预畸变处理频域验证基于带通滤波器****************************************
%  传递函数为： C =     s
%                  ————————————
%                   s² + wc²

%！ **不进行双线性变换的预畸变处理**
% 连续系统的带通中心频率为 Wc
% 但离散系统的中心频率发生的畸变，
wc = 100; % 中心频率Hz
T = 0.002;
Wc = 2 * pi * wc; % 中心频率rad/s
Ca = tf([1,0],[1 0 Wc^2]);  
Cd = c2d(Ca, T,'tustin');
figure(7)
bode(Ca, Cd);
legend('Ca','Cd')

%! **进行预畸变处理**
% 使得预畸变后的连续系统通过双线性变换后的中心频率在我们期望的Wc处
Wc_deal = 2/T * tan(Wc*T/2);
Ca_deal = tf([1, 0], [1 0 Wc_deal^2]);
Cd_deal = c2d(Ca_deal, T,'tustin');
figure(8)
bode(Ca_deal, Cd_deal);
legend('Ca_deal', 'Cd_deal');
% **********************************************************************************************************

参考文章

当然我是基于我的理解去表达的，读者可能或多或少存在这样那样疑问，我筛选出我看过质量比较好的文章附上
[1].模拟滤波器转化为数字滤波器
[2]. 双线性变换
[3]. 常用的传递函数离散化方法，差分，双线性，预矫正双线性，matlab源码
[4]. 数字信号处理——IIR滤波器设计
[5]. 使用Matlab设计数字滤波器，从原理到代码