Homework 1 AerE 432 Fall 2016 Due 9/7(W) SOLUTION

Problem 1 (20pts) The purpose of this problem is to investigate the strategy that is used in Matlab to sample a continuous-time impulse response associated with a 1st order system. Consider a first order system . It should be clear that this system has a static gain which is and a break frequency . When using the Matlab command impulse,m for this system transfer function, you can also obtain the sampling interval that was used. Find this sampling interval for each of the following values of , in order to determine how Matlab selects the sampling frequency as a function of the break frequency.

Solution: [See Appendix code.]

From the impulse command, we find that the sampling frequencies are: ws=13.6 136 1364 r/s. {In the 2015 version they were {ws = 5.9294 59.2942 592.9416. So Matlab decided to double them. Why?}

Hence, for a first order system Matlab samples the impulse response at ~ 136 times the break frequency. Hence, the Nyquist frequency is . At this frequency the FRF magnitude is

.

Remark. It can be observed that the Bode plots go to -40dB. How is it computed?


Problem 2 (40pts) The purpose of this problem is to continue the investigation for the case of a 2nd order under damped system

This system has unity static gain, an undamped natural frequency and a damping ratio .

(a)(20pts)) Assuming this system has a sufficiently small damping ratio, then the peak in the FRF occurs at the undamped natural frequency. Find the magnitude of this peak. Also, find the half-power frequency. Both of your answers should be in terms of .

Solution:

(i) . Hence, the magnification at resonance is or

for small .

(ii) By definition, a half-power frequency is a frequency ω such that [Notice that power is interpreted as the square of amplitude.] Hence,

. This relation gives , which can be written as , which is quadratic in the variable . Hence, from the quadratic formula, we have

We can simplify this as follows:

.

Since , we have . And so: the hlaf-power bandwidth has

and .

(b)(20pts)) For each of the values perform the procedure used in 1. above to obtain the sampling frequency. Then attempt to determine how Matlab chooses this. Is it a multiple of the half-power frequency? The natural frequency?

ζ / 0.1 / 0.01 / 0.001 / 0.0001
ωs [rad/sec] / 20 / 20 / 20 / 2.73

Solution: The results are shown in the table at right.

For all ζ-values but 0.0001 Matlab samples at ~20r/s.

Bode plots show that the Nyquist frequency of ~10r/s

corresponds to the -40dB level. For the 0.0001 value

the sampling frequency is , and so the

Nyquist frequency is . At this frequency

the magnitude is ~1.37dB; not -40dB. It is, however,

~72.6dB below the peak magnitude. Overlaying

the impulse response for and ,

the sampled system has a larger magnitude than the analog system;

as evidenced in the zoomed plot below.

Conclusion: For very low damping Matlab seems to change from the -40dB re: static gain rule to a rule that relates to the peak magnification. It is interesting that it uses the -37dB rule for first order systems and a 40dB rul for most second order underdamped ones.

Problem 3 (40pts) The purpose of this problem is to investigate how sampling influences the FRF magnitude of a continuous-time system.

(a)(20pts) For a first order system with static gain one, and break frequency 1 r/s plot the theoretical FRF magnitude on log-frequency scale extending from 0.01 to 100. Then overlay the FRF magnitude associated with the sampling frequency used by Matlab to arrive at the system impulse response.

There are a number of ways to compute the FRF of the sampled system. Here, I will obtain it via the sampled system transfer function. Since where

. Hence, .

Remark: If we ‘normalize’ ω via [rad], then we have . Many books take this approach. I personally prefer to not normalize frequency.

(b)(20pts) Repeat (a) for a second order under damped system with natural frequency one, and for the damping ratio values given in 2(b).

Solution: From the table inside the front cover of the book, we have

Transform pair: / / /

To use this, write

Then: where

Discuss how sampling “distorts” or “aliases” the true FRF magnitude (dB) information.

Appendix Matlab Code

%Program Name: hw1.m

%PROBLEM 1

wb = [0.1 1.0 10];

dt1=zeros(1,3);

ws1 = zeros(1,3);

for k = 1:3

H = tf(1,[1/wb(k) 1]);

[h,t] = impulse(H); %Including ‘t’ allows you to compute dt.

dt1(k) = t(2) - t(1);

ws1(k) = 2*pi/dt1(k);

end

ws1

pause

%======

% PROBLEM 2

dt2=zeros(1,4);

ws2 = zeros(1,4);

z = [0.1 0.01 0.001 0.0001]; %Damping Ratios

for k = 1:4

H = tf(1,[1 2*z(k) 1]);

[h,t] = impulse(H);

dt2(k) = t(2) - t(1); %Matlab sampling interval for z(k)

ws2(k) = 2*pi/dt2(k);

end

[z;ws2]

pause

%======

% PROBLEM 3

% PART (a):

% ANALOG system FRF:

figure(8)

w=0.01:0.01:100;

H = (1+i*w).^-1; %Analog system FRF

mHdB=20*log10(abs(H));

thH=atan2(imag(H),real(H))*180/pi; %Phase in degrees

subplot(2,1,1),semilogx(w,mHdB)

xlabel('Frequency (rad/sec)')

ylabel('dB')

grid

hold on

subplot(2,1,2),semilogx(w,thH)

xlabel('Frequency (rad/sec)')

ylabel('Phase (degrees)')

grid

hold on

%------

% SAMPLED system:

wsamp =ws1(2); %From Problem 1

wN =wsamp/2; %Nyquist frequency

del = 2*pi/wsamp;

a = exp(-del);

ww=.01:.01:wN;

z1=exp(-1i*ww*del);

HH = del *(1-a*z1).^-1;

mdB = 20*log10(abs(HH));

subplot(2,1,1),semilogx(ww,mdB,'r')

grid

title('1ST Order Analog (blue) & Sampled (red) FRF Magnitudes')

th=atan2(imag(HH),real(HH))*180/pi;

subplot(2,1,2),semilogx(ww,th,'r')

grid

title('1ST Order Analog (blue) & Sampled (red) FRF Phases')

pause

%======

% PART (b)

for k = 1:4

sys = tf(1 , [1 2*z(k) 1]);

figure(10+k)

%ANALOG System:

w=0.01:0.1:100;

H = ((1i*w).^2 + 2*z(k)*(1i*w) + 1).^-1; %Analog system FRF

mHdB = 20*log10(abs(H));

subplot(2,1,1),semilogx(w,mHdB)

grid

hold on

thH=atan2(imag(H),real(H))*180/pi;

subplot(2,1,2),semilogx(w,thH)

grid

hold on

%SAMPLED System:

wsamp =ws2(k);

wN=wsamp/2;

ww=.01:.001:wN;

del = 2*pi/wsamp;

a = exp(-del*z(k));

wd = (1-z(k)^2)^0.5;

c = del/(1-z(k)^2);

ww = 0.01:0.001:wN;

z1 = exp(-1i*ww*del);

num = c*a*sin(wd*del)*z1;

den=1 - 2*a*cos(wd*del)*z1 + a^2*z1.^2;

HH=num./den; %Sampled system FRF

mdB = 20*log10(abs(HH));

subplot(2,1,1),semilogx(ww,mdB,'r')

xlabel('Frequency (rad/sec)')

ylabel('dB')

title(['Analog (blue) & Sampled (red) FRF Magnitudes for',num2str(z(k)),' '])

th=atan2(imag(HH),real(HH))*180/pi;

subplot(2,1,2),semilogx(ww,th,'r')

xlabel('Frequency (rad/sec)')

ylabel('dB')/

title(['Analog (blue) & Sampled (red) FRF Phases for',num2str(z(k)),' '])

end