# Sinusoidal Frequency Response of a Transfer Function

42 views (last 30 days)
Mustafa Furkan on 15 Jan 2023
Edited: Paul on 15 Jan 2023
I have a transfer function as below. I want to plot the (FRF) response of a function like y(t)=A.sin(w.t) for this transfer function. The A and w values are constants, any value can be substituted. I think the response will be steady-state.
A = [0,1,0,0;-64,0,64,0;0,0,0,1;64,0,-128,0];
B = [0;0;0;64];
C = [1,0,0,0]; D = 0;
sys1 = ss(A,B,C,D);
sys2 = tf(sys1)
sys2 = 4096 -------------------------------------------------- s^4 - 1.554e-15 s^3 + 192 s^2 - 6.977e-14 s + 4096 Continuous-time transfer function.
##### 2 CommentsShowHide 1 older comment
Mustafa Furkan on 15 Jan 2023
Edited: Mustafa Furkan on 15 Jan 2023
@Paul I'm doing analysis on a system like in the image Here y represents the input. I want to plot the response of the system based on this input value. As I mentioned in the question, I want to find a solution according to an equation like y(t)=A.sin(wt). Since it is coupling, steady state will be the right choice in terms of obtaining cross spectral density. That's why I specifically mentioned the steady state in the question.

Star Strider on 15 Jan 2023
Edited: Star Strider on 15 Jan 2023
If I understand your Question correctly, the bode, bodeplot, nichols, freqresp and lsim functions will do what you want, or for a specific frequency and amplitude, the lsim funciton is also an option —
A = [0,1,0,0;-64,0,64,0;0,0,0,1;64,0,-128,0];
B = [0;0;0;64];
C = [1,0,0,0]; D = 0;
sys1 = ss(A,B,C,D);
sys2 = tf(sys1)
sys2 = 4096 -------------------------------------------------- s^4 - 1.554e-15 s^3 + 192 s^2 - 6.977e-14 s + 4096 Continuous-time transfer function.
figure
stepplot(sys1) s = stepinfo(sys1)
s = struct with fields:
RiseTime: NaN TransientTime: NaN SettlingTime: NaN SettlingMin: NaN SettlingMax: NaN Overshoot: NaN Undershoot: NaN Peak: Inf PeakTime: Inf
t = linspace(0, 100, 1001);
Ts = (t(2)-t(1))
Ts = 0.1000
Fs = 1/Ts
Fs = 10
ufcn = @(t,A,omega) A.*sin(omega.*t);
u = ufcn(t, 2.5, 150);
y = lsim(sys1, u, t);
figure
plot(t, y)
grid
xlabel('Time (s)')
ylabel('Amplitude (units)') figure
plot(t, y)
grid
xlim([0 10])
xlabel('Time (s)')
ylabel('Amplitude (units)') The ‘sys1’ and ‘sys2’ functions are obviously the same in this application. I just used them to illustrate the two plotting functions.
Further information is available in What is Frequency Response?
EDIT — (15 Jan 2023 at 21:46)
Eliminated everything except the lsim code, added stepplot plot and stepinfo call.
I doubt that a steady state is possible with this system.
.
Star Strider on 15 Jan 2023
My pleasure!
Yes —
A = [0,1,0,0;-64,0,64,0;0,0,0,1;64,0,-128,0];
B = [0;0;0;64];
C = [1,0,0,0]; D = 0;
sys1 = ss(A,B,C,D);
sys2 = tf(sys1)
sys2 = 4096 -------------------------------------------------- s^4 - 1.554e-15 s^3 + 192 s^2 - 6.977e-14 s + 4096 Continuous-time transfer function.
omegav = linspace(0, 300, 5000);
figure
bode(sys1, omegav)
grid [~,~,omega] = bode(sys1, omegav);
omegalimits = [min(omega) max(omega)]
omegalimits = 1×2
0 300
The ‘omegav’ vector defines the radian frequency vector that this bode call uses to calculate the frequency response. The frequency axis scale is logarithmic, however it goes from 0 to 300 rad/s, as the ‘omega’ result demonstrates.
.

Paul on 15 Jan 2023
If the system were bounded-input-bounded-output (BIBO) stable, then the steady state output in response to input y(t) = A*sin(w*t) would be zss(t) = M*A*sin(wt + phi), where M and phi are determined by the magnitude and phase of the system transfer function evaluated at s = 1j*w.
However, this system has no damping and is not BIBO stable, as indicated by the eigenvalues of the A-matrix all on the imaginary axis
A = [0,1,0,0;-64,0,64,0;0,0,0,1;64,0,-128,0];
format long e
e = eig(A)
e =
0.000000000000000e+00 + 1.294427190999917e+01i 0.000000000000000e+00 - 1.294427190999917e+01i 0.000000000000000e+00 + 4.944271909999160e+00i 0.000000000000000e+00 - 4.944271909999160e+00i
So, in this case we have to consider two cases:
1. w not equal to 12.94427.... and w not equal to 4.94427.... In this case, the output of the system will be the sum of zss as defined above and and two other sine waves at frequencies 12.944... an 4.94427.... The amplitude and phase of these two other sine waves would have to be determined based on A and w.
2. w = 12.9442.... or w = 4.94427.... In this case the output will include a term like C*t*sin(w*t), i.e., the output will grow indefinitely
Finding closed form expressions in either case is feasible, but could be painful.
In either case, you could use lsim to approximate the output via simulation; make sure the the time vector has a time separation small relative to the larger of 12.944 or w, mabye something like dt < 1/10/max(w,12.9444)
B = [0;0;0;64];
C = [1,0,0,0]; D = 0;
sys1 = ss(A,B,C,D);
% Case 1 example
w = 8; a = 1;
dt = 1/12.94/10;
t = 0:dt:10;
u = a*sin(w*t);
y = lsim(sys1,u,t);
plot(t,y) As expected, the output is the sum of three sinusoids at the frequencies determined by w and the imaginary part of e.
Y = fft(y);
w = (0:(numel(Y)-1))/numel(Y)*2*pi/dt;
figure
plot(w,abs(Y))
xlim([0 20]) % Case 2 example
w = abs(imag(e(1))); a = 1;
dt = 1/12.94/10;
t = 0:dt:10;
u = a*sin(w*t);
y = lsim(sys1,u,t);
figure
plot(t,y),grid This output is dominated by the growing sine wave at 12.444 rad/sec, but it does have a 4.944 rad/sec component
Y = fft(y);
w = (0:(numel(Y)-1))/numel(Y)*2*pi/dt;
figure
plot(w,abs(Y))
xlim([0 20]) This analysis would be much simpler if you added a dashpot (or damper) in parallel with two springs.
##### 2 CommentsShowHide 1 older comment
Paul on 15 Jan 2023
Edited: Paul on 15 Jan 2023
Use bode or bodeplot to compute/plot the frequency response from the state space model.

### Categories

Find more on Time and Frequency Domain Analysis in Help Center and File Exchange

R2022b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!