%図11.2
%
clear all; close all;

%1次遅れ系のパラメータ
T = 1; %T=1
K = 1; %K=1

%伝達関数の分子・分母多項式を与える
num = [ 0 K ]; %分子多項式
den = [ T 1 ]; %分母多項式

%伝達関数表現を与える
sys = tf( num, den ); %伝達関数表現

%図11.2(a)
%正弦波A sin ωtのパラメータを与える
A = 1; %正弦波の振幅
omega1 = 0.01; %ω=0.01の場合の正弦波の角周波数

%時間変数の定義
t1 = 0:0.1:1000; %0から1000まで0.1刻み

%正弦波(入力)の計算
u1 = A*sin( omega1*t1 ); %ω=0.01の正弦波を与える

%応答の初期値を与える
x0 = [ 0 ]; %x(0)=0と与える

%周波数応答(出力)を求める
[ y1, t1, x1 ] = lsim( sys, u1, t1, x0 ); %周波数応答

%図11.2(a)のプロット

figure(1) %図のウィンドウを開く

plot(t1,u1,'-b',t1,y1,'-r'); %周波数応答をプロット
xlim([0,1000]) %横軸(時間軸)の範囲の指定)
ylim([-1.0 1.0]) %縦軸の範囲の指定
xticks([0 200 400 600 800 1000]) %横軸の目盛りの値の設定
yticks([-1.0 -0.5 0.0 0.5 1.0]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('input and output') %縦軸のラベル表示
legend('input','output'); %凡例の表示

%図11.2(b)
%正弦波A sin ωtのパラメータを与える
A = 1; %正弦波の振幅
omega2 = 10; %ω=10の場合の正弦波の角周波数

%時間変数の定義
t2 = 0:0.01:10; %0から10まで0.01刻み

%正弦波(入力)の計算
u2 = A*sin( omega2*t2 ); %ω=10.0の正弦波を与える

%応答の初期値を与える
x0 = [ 0 ]; %x(0)=0と与える

%周波数応答(出力)を求める
[ y2, t2, x2 ] = lsim( sys, u2, t2, x0 ); %周波数応答を求める

%図11.2(b)のプロット
figure(2) %図のウィンドウを開く
plot(t2,u2,'-b',t2,y2,'-r'); %周波数応答をプロット
xlim([0,10]) %横軸(時間軸)の範囲の指定
ylim([-1.0 1.0]) %縦軸の範囲の指定
xticks([0 2 4 6 8 10]) %横軸の目盛りの値の設定
yticks([-1.0 -0.5 0.0 0.5 1.0]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベルの表示
ylabel('input and output') %縦軸のラベル表示
legend('input','output'); %凡例の表示

%図11.2(c)のプロット

figure(3) %図のウィンドウを開く

plot(t2,u2,'-b',t2,y2,'-r'); %周波数応答をプロット
xlim([5,10]) %横軸(時間軸)の範囲の指定
ylim([-1.0 1.0]) %縦軸の範囲の指定
xticks([5 6 7 8 9 10]) %横軸の目盛りの値の設定
yticks([-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0]) %縦軸の目盛りの値の設定
grid; %罫線を表示
xlabel('time  t [s]'); %横軸のラベル表示
ylabel('input and output') %縦軸のラベル表示
legend('input','output'); %凡例の表示