python 频率响应 反求 传递函数

一、理论公式

        频率响应函数的表达式：

        对应的z变换的多项表达式：

        Z变换的零极点表达形式：

        Z变换的二阶因子级联形式：

二、滤波函数filter

1. filter函数，仅可以用于零状态响应系统。

        y=filter(b,a,x) ；       

        %b为z变换多项表达式公式中[b0,b1...bM]的矩阵

        %a为z变换多项表达式公式中[a0,a1...aN]的矩阵

        %x为输入的原始信号

移动平均滤波的matlab程序示例：

%求解零状态差分方程函数--LTI系统
% y=filter(b,a,x);


N=201;
n=rand(1,N)-0.5;%噪声信号
k=0:N-1;%
x=2*k.*(0.9.^k)+2.0*cos(0.02*pi*k)+5.0;%输入信号
y=x+n;%包含噪声的信号
plot(k,n,'g--',k,x,'b--',k,y,'r-');%显示三个函数图像，
xlabel('Time index k');
legend('n[k]','x[k]','y[k]');%添加图例


M=10;
b=ones(M,1)/M;%移动平均滤波，10个点相加，再除10，算出平均值
a=[1];
filter_y=filter(b,a,y);
plot(k,x,'b-',k,y,'g--',k,filter_y,'r-');%显示三个函数的波形
xlabel('Time index k');
legend('x[k]','y[k]','filter_y[k]');%添加图例

得出的波形如下图：

三、频率响应函数

 H=freqz（b,a,w）；

%b为z变换多项表达式公式中[b0,b1...bM]的矩阵

        %a为z变换多项表达式公式中[a0,a1...aN]的矩阵

        %w为输入的角频率

matlab程序示例

%freqz用以分析离散系统的频率响应

b=[1]; %分子矩阵
a1=[1,-0.9]; a2=[1,0.9];%分母矩阵
w=linspace(0,pi,512);%在0-π范围内等间隔分512份
h1=freqz(b,a1,w);%计算频率响应
h2=freqz(b,a2,w);%计算频率响应
% plot(w,abs(h1),w,abs(h2),':');
plot(w/pi,abs(h1),w/pi,abs(h2),':');
legend('a=0.9','a=-0.9');

 

可以看出在分母矩阵为[1,-0.9]时，系统为低通滤波器；分母矩阵为[1,0.9]时，系统为高通滤波器。

 四、频率响应函数不同形式的转换

b=[1,4];%Z变换的分子矩阵
a=[1,0.1,-0.2];%Z变换的分母矩阵
z=zeros(1,2);%2*1的矩阵，值为0
% [z,p,K]=tf2zp(b,a);%多项式转零极点表达式
[z,p,K]=tf2zpk(b,a);%多项式转零极点表达式
sos=zp2sos(z,p,K);%零极点转二阶因子级联形式

多项式表达式为

经过matlab计算，得到其零极点表达式为

二阶因子级联形式：

         （只有1级）

需要注意的是：tf2zp计算出来的结果中只有一个零点，而 tf2zpk计算出来的结果中有两个零点（多出来一个位置在0处的零点）

        tf2zpk是用于离散系统传递函数，基于z变换；而tf2zp是用于连续系统传递函数，基于拉普拉斯变换。依据二者的函数公式，计算该系统函数的零极点出来，确实会存在零点的细微差别。此处不做讨论，注意即可。