【模拟滤波与数字滤波区别】--双线性变换的预畸变处理理解及验证；matlab中freqz、freqs、loglog、bode函数绘图的区别
%% 滤波器练习
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。
% 在设计离散滤波器的时候进行了预畸处理，所以其离散系统bode图截止频率在截止频率30Hz处，同时对应的连续系统Wc不在30Hz处
% 角频率转换关系：Ω = 2/T *tan(W*T/2)。所以当W*T/2足够小时，Ω=W，为线性转换关系，基于bilinear转换的离散和连续系统的失真比
% 较小；另一个角度可以看出T过大或者W过大，都会加剧失真
% 得到的频率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');
% **********************************************************************************************************