扫频法求开环传递函数_传递函数扫频_MATLAB程序_扫频法求传递函数

% Transfer function identification with frequency test 扫频测试法辨识传递函数，再计算得到传递函数； clear all; close all; load idenfile; %读取扫频数据 a=25;b=133;c=10; %传递函数参数值 sys=tf(b,[1,a,c]); %构造传递函数sys 这个是原来的传递函数； ts=0.001; %离散化的时间； Am=0.5; %幅值 kk=0; %变量kk=0 q=1; %变量q=1 for F=1:0.5:10 %循环，实际是品绿循环，从1开始，每隔0.5 到10 这个是扫频的频率值 kk=kk+1; % 1,1.5,2,2,5......10 共有19个； kk = 19; FF(kk)=F; % FF() = F; for i=10001:1:15000 % 求取数据是20000个数据的10001到15000； fai(1,i-10000) = sin(2*pi*F*i*ts); % 求取“φ”1 的实部 fai 1 是实部 fai(2,i-10000) = cos(2*pi*F*i*ts); % 求取“φ”2 的虚部 fai 2 是虚部； end Fai=fai'; % 矩阵共轭转置，fai.' 非共轭转置模式，fai'共轭转置，当为实矩阵时候，相同 % 单纯地共轭用：conj() % 单纯的转置用：transpose() % 转换完成后，Fai 就变成了 2x5000的矩阵；每一行sin,cos fai_in(kk)=0; % 定义输入的角“φ”， Y_out=y(q,10001:1:15000)'; % 定义一个Y_out 矩阵，再非共轭转置，y()是一个19x20000的矩阵；q会累加 cout=inv(Fai'*Fai)*Fai'*Y_out; % inv：Y=inv(X)返回方阵X的逆矩阵 fai_out(kk)=atan(cout(2)/cout(1)); % Phase Frequency(Deg.) 相频（度） if fai_out(kk)>0 % 如果计算的值 > 0 fai_out(kk)=fai_out(kk)-pi; % 需要减去pi; end Af(kk)=sqrt(cout(1)^2+cout(2)^2); % Magnitude Frequency(dB) 得到幅值； mag_e(kk)=20*log10(Af(kk)/Am); % in dB. 幅值比； ph_e(kk)=(fai_out(kk)-fai_in(kk))*180/pi; % in Deg. 相位差； if ph_e(kk)>0 ph_e(kk)=ph_e(kk)-360; end q=q+1; end %以上是根据扫频数据计算得到幅值比和相位差的值，一共19个，从1，间隔0.5 到10 %====================================================== FFF=FF'; % 再讲FF 变换转置 ，在这里我将FF 改变为FFF %%%%%%%%%%%%%%% Closed system modelling %%%%%%%%%%%%%%%%%% mag_e1=Af'/Am; %From dB.to ratio ph_e1=fai_out'-fai_in'; %From Deg. to rad hp=mag_e1.*(cos(ph_e1)+j*sin(ph_e1)) ; %Practical frequency response vector na=2; % Second order transfer function 设置为2阶的传递函数； nb=0; w=2*pi*FFF; % in rad./s 原来是FF 我给改成了FFF % bb and aa gives real numerator and denominator of transfer function bb和aa给出了传递函数的实分子分母 [bb,aa]=invfreqs(hp,w,nb,na); % w(in rad./s) contains the frequency values 从扫频数据得到传递函数； bb % 显示bb 值 aa %显示 aa 值 G=tf(bb,aa) % Transfer function fitting 建立一个传递函数G %---------------------------------------------------------------------------------------------------------------------------------------- %sys=tf(b,[1,a,c]); % tf（）是建立一个系统的传递函数，根据上面给定的传递函数值，建立sys变量的传递函数 % tf() 是传递函数，即被控对象函数G(); dsys=c2d(G,ts,'z'); % c2d（）从s域转换为z域，离散化，步长为ts=0.001; % 把控制函数离散化(离散步长为0.001)后，取Z变换n阶定常离散系统查分方程； % dsys 即是 Y(z)/U(z) !!! [num,den]=tfdata(dsys,'v'); % 加上'v'，可以让输出的bai值由元胞数组du改为数组直接输出： % 离散化后提取查分方程的分子分母； num den %---------------------------------------------------------------------------------------------------------------------------------------- hf=freqs(bb,aa,w); % Fited frequency response vector % hf = freqs(b, a, w) 根据系数向量计算返回模拟滤波器的复频域响应。 % freqs(b, a, f) 挑选f个频率点来计算频率响应h % Transfer function verify: Getting magnitude and phase of Bode %------------------------------------------------------------------------------------ %下面这几句话和以前的一样！！！ sysmag=abs(hf); % ratio. sysmag1=20*log10(sysmag); % From ratio to dB sysph=angle(hf); % Rad. sysph1=sysph*180/pi; % From Rad.to Deg. %------------------------------------------------------------------------------------ % Compare practical Bode and identified Bode figure(1); subplot(2,1,1); semilogx(w,mag_e,'r',w,sysmag1,'b');grid on; xlabel('rad./s');ylabel('Mag.(dB.)'); subplot(2,1,2); semilogx(w,ph_e,'r',w,sysph1,'b');grid on; xlabel('rad./s');ylabel('Phase(Deg.)'); figure(2); subplot(2,1,1); magError=sysmag1-mag_e'; plot(w,magError,'r'); xlabel('rad./s');ylabel('Mag.(dB.)'); subplot(2,1,2); phError=sysph1-ph_e'; plot(w,phError,'r'); xlabel('rad./s');ylabel('Phase(Deg.)'); figure(3); bode(sys,'r',G,'b');