画单自由度系统传递函数(实频，虚频，幅频,相位，导纳)

激励采用脉冲激励，脉冲激励为Hanning函数

function f = hanning_imp(t, Tc, A)
f = zeros(size(t));
f(t < Tc) = A / 2 * (1 - cos(2*pi * t(t < Tc) / Tc));
end

下面是计算单自由度系统响应的Matlab程序，计算传递函数，画实频、虚频、幅频、相位、导纳图

m = 100;
k = 1000;
c = 100;

num = 1;
den = [m c k];
sys = tf(num, den);
dt = 0.00001;
fs = 1/dt;   %采样频率(Hz) 100Hz 实际并不需要这么高的采样频率，但是如果采样时间太小，hanning脉冲不完整
% 为了得到准确的响应dt一定要小，否则做出的相位可能不对

t = 0:dt:200;
Tc = 0.001;
A = 10;
u = hanning_imp(t, Tc, A);
y = lsim(sys, u, t);
% y = impulse(sys, t);
y = y';

N = length(u);
fy = fft(y);
fu = fft(u);
ft = fy ./ fu;
f = (0:N-1) * fs ./ N;
ft_r = real(ft);
ft_i = imag(ft);

part = (f < 30);

figure('name', '实频');
plot(f(part), ft_r(part));
figure('name', '虚频');
plot(f(part), ft_i(part));
figure('name', '幅值');
plot(f(part), abs(ft(part)));
figure('name', '相位');
plot(f(part), phase(ft(part)) * 180 / pi);
figure('name', '导纳圆');
plot(ft_r(part), ft_i(part), '.');
axis equal;
hold on;
xk = [-0.0001165
0.0002857
8.297e-5
0.0004869
0.000682];
xk1 = [-0.0008017
-0.0003082
0.0002857
0.001194
0.001451];
yk1 = [-0.002817
-0.003095
-0.003185
-0.002869
-0.002633];

yk = [-0.003149
-0.003185
-0.003179
-0.003165
-0.003121];
k = 5;
A = [
sum(xk.^2)  sum(xk .* yk)   sum(xk)
sum(xk.*yk) sum(yk.^2)   sum(yk)
sum(xk)  sum(yk) k];
B = -[
sum(xk.^3 + xk.* yk.^2)
sum(xk.^2.*yk + yk.^3)
sum(xk.^2+yk.^2)];
rlt = A\B;
x0 = rlt(1)*(-0.5);
y0 = rlt(2)*(-0.5);
r = sqrt(rlt(1)^2/4+rlt(2)^2/4-rlt(3));
fai = 0:0.01:2*pi;
x = x0 + r * cos(fai);
y = y0 + r * sin(fai);
plot(x,y, 'r');

（1） 导纳图（还有拟合的圆，采用的拟合数据没有更新，不是这个导纳圆的）



（2）幅频图



（3）相频图



（4）实频图



（5）虚频图