Matlab+homework+2

Part A (In-class room exercises)

1. Using plot to display the following voltage with appropriate line type, title and labels.

Also present the graph with suitable ranges of axis.

v(t)220sin2ft

where f=50 Hz.

(Hint: time interval must be small enough, e.g. 0.001 seconds. Therefore, t=0:0.001:1 is appropriate);

t=0:0.002:1;

v=220*sin(2*pi*50*t); plot(t,v,'r','LineWidth',1); xlabel('time'); ylabel('value'); title('v(t)--t');

axis([0 1 -250 250]);

2. Using plot to display the following voltage with appropriate line type, title and labels.

Also present the graph with suitable ranges of axis.

v(t)220e5tsin2ft

where f=50 Hz.

In addition, on the same graph, draw the envelope of the oscillation and add legends. t=0:0.001:0.1;

v=220*sin(2*pi*50*t);

v1=220*exp(-5*t).*sin(2*pi*50*t); plot(t,v,'r','LineWidth',1); hold on;

plot(t,v1,'b','LineWidth',2); xlabel('time'); ylabel('value'); title('v(t)--t');

axis([0 0.1 -250 250]); hold off; grid on;

3. Use subplot, draw a 2 by 2 array of plots for the following functions:

v1cos10t,v2e5tcos10tv3e10tcos10tv4e20tcos10t

Apply appropriate line type, title, labels and axis ranges for the graphs. t=0:.01:1

v1=cos(10*pi*t);

v2=exp(-5*t).*cos(10*pi*t); v3=exp(-10*t).*cos(10*pi*t); v4=exp(-20*t).*cos(10*pi*t); subplot(2,2,1); plot(t,v1); xlabel('time'); ylabel('value'); title('v--t');

axis([0 1 -1 1]); grid on;

subplot(2,2,2); plot(t,v2); xlabel('time'); ylabel('value'); title('v--t');

axis([0 1 -1 1]); grid on;

subplot(2,2,3); plot(t,v3); xlabel('time'); ylabel('value'); title('v--t');

axis([0 1 -1 1]); grid on;

subplot(2,2,4); plot(t,v4); xlabel('time'); ylabel('value'); title('v--t');

axis([0 1 -1 1]); grid on;

4. Use plot3 to plot 2 spiral curves like below with appropriate line width and colour.

t=0:0.01:10;

x=0.01*sin(2*pi.*t).*exp(0.6*t); y=0.01*cos(2*pi.*t).*exp(0.6*t); z=0.2*t;

plot3(x,y,z,'g','LineWidth',3); axis([-4 4 -4 4 0 3]);

grid on;

32.521.510.50420-2-4-4-2204

t=0:0.01:10;

x=sin(2*pi.*t).*log(6*t); y=cos(2*pi.*t).*log(6*t); z=0.2*t;

plot3(x,y,z,'g','LineWidth',3); axis([-4 4 -4 4 0 3]); grid on;

32.521.510.50420-2-4-4-2204

5. Display the surface using mesh and contour with a suitable resolution:

0.005((x50)z(x,y)e

2(y50)2)

[x,y]=meshgrid(0:2:100,0:2:100); z=exp(-0.005*((x-50).^2+(y-50).^2)); figure(1); mesh(x,y,z), xlabel('x'); ylabel('y'); grid; figure(2); contour(x,y,z); xlabel('x');

ylabel('y'); grid;

10.80.60.40.2010010050020040x6080y

100908070605040302010001020304050x60708090100y

6. Load 2 of your photos into Matlab WorkSpace using imread. a) Change brightness

locally or globally. b) Overlap them to produce a new photo. c) write a new photo into a file using imwrite. (Note: lower version of Matlab such as 6.5 is not allowed to do some direct image operations.) A):

A=imread('Desert.jpg'); figure(1),imshow(A); size(A);

figure(2),histeq(A(:,:,1)) imshow(A+100);

imwrite(A+100,'ddd1.jpg');

B):

A=imread('Desert.jpg'); B=imread('Penguins.jpg'); [M,N,L]=size(A);

C(1:M,1:N,1:L)=B(1:M,1:N,1:L); figure(1);imshow(A); figure(2);imshow(C); figure(3);imshow(A+C);

C):

A=imread('Tulip.jpg'); imwrite(A,'ttt.jpg');

7. For linear simultaneous equations

a11x1a21x1...a12x21...a22x2......a1NxNa2NxN...b1b2...

aMx1aM2x2...aMNxNbMthe equation coefficients: A=[ 1 -1 4 3 -5 4 5 -6 0 7 -8 9

-1 3 -2 6]; (M=N)

a) Find the determinant of A, >> A=[1 -1 4 3; -5 4 5 -6; 0 7 -8 9 ; -1 3 -2 6]; >> det(A)

ans =

-765.0000

b) Find the inverse of A and check for matrix singularity, >> inv(A)

ans =

0.3333 -0.0000 0.3333 -0.6667 0.2000 0.1176 0.2706 -0.3882 0.2000 0.0588 0.0353 -0.0941 0.0222 -0.0392 -0.0680 0.2183

c) If B=[ 5; 1; -2; 3], find the unknown x in the equation. >> B=[5;1;-2;3]; >> x=inv(A'*A)*A'*B x =

-1.0000

-0.5882 0.7059 0.8627

8. For linear simultaneous equations,

a11x1a21x1...a12x21...a22x2......a1NxNa2NxN...b1b2...

aMx1aM2x2...aMNxNbMM>N, the equation coefficients: A=[ 1 -1 4 3 -5 4 5 -6 0 7 -8 9 -1 3 -2 6 1 -2 5 3

0 4 7 -3];

a) Find the determinant of A’*A, >> A=[1 -1 4 3; -5 4 5 -6; 0 7 -8 9; -1 3 -2 6; 1 -2 5 3; 0 4 7 -3]; >> det(A'*A)

ans =

2.1051e+007

b) Find the inverse of A’*A and check for matrix singularity, >> inv(A'*A)

ans =

0.0884 0.0321 -0.0013 -0.0219 0.0321 0.0231 0.0002 -0.0099 -0.0013 0.0002 0.0085 0.0053 -0.0219 -0.0099 0.0053 0.0144

c) If B=[ 5; 1; -2; 3; 4; 0], find the solution of the equation. >> B=[5;1;-2;3;4;0]; >> x=inv(A'*A)*A'*B

x =

-0.03 -0.48 0.5755 0.7083

9. Data of 10 records are shown below

y=[3.5 4.3 3.7 5.4 6.6 7.3 8.7 8.8 9.4 9.0 10.0 12.0 11.3 9.9 13.3],

Use polyfit with different orders (from 1 to 3) of polynomials to find a curve of best fit. Check the total distance between the fitted curve z and records defined by

2szyiiix=0:14;

1/2.

y=[3.5 4.3 3.7 5.4 6.6 7.3 8.7 8.8 9.4 9.0 10.0 12.0 11.3 9.9

13.3];

plot(x,y,'x'); hold on;

p=polyfit(x,y,3); z=polyval(p,x);

plot(x ,z ,'r', 'LineWidth', 2); hold off; grid on;

1412108202468101214

>> s=sqrt(sum((z-y).^2)) s =

3.0398

10. Create a set of 20 points from a curve by Matlab code:

x=1:20;y=2*exp(-0.3*(x-5).^2)+0.7*exp(-0.2*(x-12).^2);

Then interpolate the curve to 60 points using ‘linear’ and ‘spline’ options, respectively. See the quality of different types of interpolation.

x=1:20;

y=2*exp(-0.3*(x-5).^2)+0.7*exp(-0.2*(x-12).^2); xi = 0:.25:20;

yi = interp1(x,y,xi,'linear'); plot(x,y,'r*'); hold on;

plot(xi,yi,'o');

plot(xi,yi,':'); hold off; grid on;

2.521.510.5002468101214161820

x=1:20;

y=2*exp(-0.3*(x-5).^2)+0.7*exp(-0.2*(x-12).^2); xi = 0:.25:20;

yi = interp1(x,y,xi,'spline'); plot(x,y,'r*'); hold on;

plot(xi,yi,'o'); plot(xi,yi,':'); hold off;

grid on;

2.521.510.50-0.502468101214161820

Part B

1. Using the plot and subplot functions create 4 plots on a 2 by 2 array of subplots,

for the function exp(-t)sin(5t) showing in each plot the function in the corresponding intervals of t i.e. (-2, -1), (-1,0), (0,1) and (1,2).

Answer:

t=-2:.01:2;

x=(exp(-t)).*sin(5.*t); %

subplot(2,2,1); plot(t,x);

axis([-2 -1 -5 5]); subplot(2,2,2);

plot(t,x);

axis([-1 0 -5 5]); subplot(2,2,3); plot(t,x);

axis([0 1 -5 5]); subplot(2,2,4); plot(t,x);

axis([1 2 -5 5]); grid on;

5500-5-25-1.5-1-5-15-0.5000-500.51-511.52

2. A three phase induction motor characteristic is given in terms of mechanical shaft output torque TM (NM Newton-meter) as a function of rotational speed ω (rad/s radian per second). This is approximated by 3 piece-wise linear equations as follows:

4.09095TM422.52010TM1200TM09090110

110120This motor is directly coupled to a load

TL5012010012032TL, which can be represented as

4120

Write 2 separate Matlab function m-files in which: a) the motor characteristic, b) the load characteristic are defined only as functions of ω. Name them motor.m and sysload.m, respectively.

Answer:

function TM=motor(w); for i=1:length(w);

if w(i)>=0 & w(i)<=90*pi; TM(i)=w(i)/(90*pi)+4.0;

elseif w(i)>=90 & w(i)<=110*pi; TM(i)=(95*w(i))/(20*pi)-422.5; elseif w(i)>=110*pi & w(i)<=120*pi; TM(i)=-10*w(i)/pi+1200; end; end; end

function TL=sysload(w); %

TL=-50*(w/(120*pi)).^3+100*(w/(120*pi)).^2+4*w/(120*pi); end

3. Write a Matlab scripgt-m file which calls your function m-files from Question 2.

And plot the motor and load characterstics on the same figure, giving suitable labelling and title. Name this script m-file systemplot.m

Answer:

w=0:1:120*pi; TM=motor(w); TL=sysload(w);

plot(w,TM,'r'); hold on;

plot(w,TL,'k'); hold off;

xlabel('rotational speed ω'); ylabel('q');

title('compare of TM & TL'); grid on;

4. Write a sript m-file which finds all the mathematicall possible points (values of ω) whereTMT L

for the range 0≤ ω≤120π (rad/s) with the characteristic given in Question 2. The system could theoretically operate at a speed ω where

TMTL Mathematically, this involves finding the roots of the equation

TMTL0

Give a name to this m-file posspoints.m . Call posspoints.m in systemplot.m and show all of the points on a system plot using the ‘o’ symbol.

Answer:

Part C

Introduction to DSP

1. Basic digital signals Unit impulse function

0n0[n]

1n0Exercise 1-1:

Display the unit impulse function with Matlab code. The code can be either typed under Matlab prompt or written into a script Matlab file, then run the file. n=-10:10; for k=1:21 x(k)=0; end; x(11)=1; stem(n, x);

axis([-10 10 0 2]);

Problem 1-1:

For the signal2[n2], write the Matlab code and copy the result figure into the following boxes. n=-10:10; 2for k=1:21 1.8 x(k)=0; 1.6end; 1.4x(9)=2; stem(n, x); 1.2axis([-10 10 0 2]); 1 0.8 0.6 0.4 0.2 0 -10-8-6-4-20246810

Unit step function

0n0u[n]

1n0

Exercise 1-2:

Display the unit step function with Matlab code:

n=-10:10; for k=1:10 x(k)=0; end;

for k=11:21 x(k)=1; end;

stem(n, x);

axis([-10 10 0 2]);

Problem 1-2:

For the signal u[n2], write the Matlab code, and display the corresponding signal into the following boxes. n=-10:10; 0for k=1:12 -0.1 x(k)=0; end; -0.2 -0.3 -0.4for k=13:21 -0.5 x(k)=-1; end; -0.6stem(n, x); -0.7axis([-10 10 0 -1]); -0.8 -0.9 -1-10-8-6-4-20246810 Sine wave

sin(nT)

where  - frequency (rad/second), T - sampling interval (seconds)

Exercise 1-3:

Let T=0.02,  =6.28. Using the following code, plot on the screen in the region 0  n 100.

n=0:100; T=0.02;

omega=6.28;

for k=1:101

x(k)=sin((k-1)*omega*T); end;

stem(n, x);

Problem 1-3:

Modify the above code to display the signal: 2sin(nT2Execute the code and plot the signal on the screen. Copy the code, then display the signal into the following boxes. n=0:100; 2T=0.02; omega=6.28; 1.5 for k=1:101 1 x(k)=2*sin((k-0.51)*omega*T+pi/2); end; 0stem(n, x); -0.5 -1 -1.5-2), where T=0.02, =6.28.

0102030405060708090100 Time shift, impulse and step responses

Exercise 1-4

For the unit impulse,

1nn0[nn0]

0nn0and the unit step

1nn0u[nn0]

0nn0write two Matlab functions and save as ‘impseq.m’ and ‘stepseq.m’ files, respectively.

% impseq.m %Unit Impulse

%[n1:n2] simple number range %n0 - time shift

function [x, n]=impseq(n0,n1,n2) n=[n1:n2]; x=[(n-n0)==0]; %end;

% stepseq.m % Unit Step

%[n1:n2] simple number range %n0 - time shift

function [x, n]= stepseq(n0, n1, n2) n=[n1:n2]; x=[(n-n0)>=0]; %end

Given the following difference equation:

y[n]0.7y[n1]0.9y[n2]x[n]

a) Calculate and plot the impulse response at n=-20,…,120; and

b) Calculate and plot the step response at n=-20,…,120, using the following code

b=[1]; a=[1,-0.7, 0.9];

x=impseq(0,-20, 120); n=[-20:120]; h=filter(b,a,x);

subplot(2,1,1);stem(n,h);

title('Impulse response');xlabel('n');ylabel('h(n)');

x=stepseq(0,-20, 120); n=[-20:120]; h=filter(b,a,x);

subplot(2,1,2);stem(n,h);

title('Step response');xlabel('n');ylabel('s(n)');

Problem 1-4 Solve the above a) and b) for the system: y[n]1.8y[n1]0.94y[n2]x[n]2x[n1]

Display the responses in the boxes. Impulse response105h(n)0-5-10-20020406080100120nStep response4030s(n)20100-2002040n6080100120 2. Convolution

In the liner time-invariant (LTI) system, the response of the system can be calculated by a convolution between the input and its unit impulse response.

Exercise 2-1 Write an Matlab function and save as ‘conv_m.m’.

%convolution function %conv_m.m

function [y, ny]=conv_m(x,nx,h,nh) nyb=nx(1)+nh(1);

nye=nx(length(x))+nh(length(h)); ny=[nyb:nye]; y=conv(x,h); %end;

Compute the convolution for x[n] and h[n] using the following codes: x=[3, 11, 7, 0, -1, 4, 2]; nx=[-3:3];

h=[2, 3, 0, -5, 2, 1]; nh=[-1:4];

[y, ny]=conv_m(x,nx,h,nh); Subplot(3,1,1); stem(nx,x); ylabel('x[n]');

axis([-6 10 -20 20]); Subplot(3,1,2); stem(nh,h); ylabel('h[n]');

axis([-6 10 -20 20]); subplot(3,1,3); stem(ny,y); xlabel('n'); ylabel('y[n]');

axis([-6 10 -50 50]);

Problem 2-1 Compute a response of a system by convolution between an input signal x1[n] and the impulse response x2[n]. Display the signals and result into the following boxes using the code modified from the above:

%********************************* %The input signal x=[1, 1, 1, 1, 1]; % ^

% The Impulse response h

=[1,0.5,0.25,0.125,0.0625,0.03125; % ^

%Define the position on axis nx=[0:4]; nh=[0:5];

%******************************* 20x[n]0-20-620-4-20246810h[n]0-20-650-4-20246810y[n]0-50-6-4-202n46810

3. Z-Transform and Frequency Response

Exercise 3-1

1. Its frequency response can be z0.81obtained by letting z =exp(j). H()

expj0.8

The MATLAB code for plotting the magnitude and phase of the frequency response:

delta=0.01;

Omega=0:delta:pi;

H= 1 ./ (exp(j .* Omega)-0.8); subplot(2,1,1);

plot(Omega, abs(H)); xlabel('0<\\Omega<\\pi'); ylabel('|H(\\Omega)|');

axis([0 pi 0 max(abs(H))]); subplot(2,1,2);

plot(Omega,atan2(imag(H),real(H))); xlabel('0<\\Omega<\\pi'); ylabel(' -\\pi < \\Phi_H <\\pi') axis([0 pi -pi pi]);

A digital LTI system has z-transform H(z)4|H()|2000.511.50<<22.530 - < H <0-1-2-300.511.50<<22.53

The LTI system has a high gain in the range of low frequencies. Therefore, it is a low-pass filter.

Problem 3-1 Generate Matlab code similar to Exercise 4-1, and sketch the magnitude

zand phase for the frequency response of a filter: H(z)2.

z1.36z0.922

15|H()|105000.511.50<<22.532 - < H <0-200.511.50<<22.53 Indicate the type of filter and comment on the magnitude and phase in different frequency bands:

Problem 3-2 Generate Matlab code similar to Exercise 4-1, and sketch the magnitude

z4H(z)2(zz0.8)(z21.36z0.922)and phase for the frequency response of a filter: 40|H()|20000.511.50<<22.532 - < H <0-200.511.50<<22.53 Indicate the type of filter and comment on the magnitude and phase in different frequency bands:

4. Spectral analysis by DFT and windowing

Digital spectral analysis using DFT – the decomposition of a signal into its frequency components using a computer or spectral-purpose hardware- is a valuable technique in many branches of engineering, applied science and information technology. An FFT algorithm is the natural choice for this work because of its fast speed. The basic assumption behind the purpose of spectral analysis is that a frequency-domain description for the signal is likely to reveal important information which is not apparent in the time–domain. Note that spectral analysis, unlike digital filtering, is primarily investigative. It is not necessarily concerned with modifying the signal. Nevertheless, the information it yields often leads to important insights or decision making.

4.1 Spectral leakage

If a component in a signal is not the harmonic of basic frequency, the discontinuity will occur between one basic sampling length and the next. The spectrum will not be represented by a single line but a peak with sidelobs.

Problem 4-1

Examine sinusoidal components in a signal by spectral analysis

232n2107.5n2422.5nx[n]0.1sin0.2sin0.15sin

102410241024

The first component gives a single spectral coefficient. The second and the third components display sidelobes around a high coefficient, which represent the spectral leakage. Using the following code, plot the signal in time and frequency domains on the screen and sketch them into the box below.

%Examine signal components subplot(2,1,1); %sampling rate 1024Hz plot(t,x); N=1024; axis([0 1 -1 1]); dt=1/N; t=0:dt:1-dt; xlabel('nT (seconds)'); ylabel('x[n]'); %****define the signal** X=fft(x); x=0.1*sin(2*pi*32 .*t ) + ... df=1; 0.2*sin(2*pi*107.5 .*t )+ ... f=0:df:N-1; 0.15*sin(2*pi* 422.5 .*t ) ; subplot(2,1,2); %*********************** plot(f,abs(X));

axis([0 N/2 0 120]); xlabel('k (Hz)'); ylabel('|X[k]|'); 10.5

x[n]0-0.5-100.10.20.30.40.50.6nT (seconds)0.70.80.91806040200020040060080010001200 4.2 Windowing

In spectral analysis, only a signal component at an exact harmonic frequency of DFT gives rise to a single and well-defined spectral line. Unfortunately practical digital signals do not behave in this way. They normally contain a wide mixture of frequencies, few of them exact harmonics. This means that the spectral leakage is generally unavoidable, and it may lead to difficulties of interpretation. It is therefore common practice to taper the original signal before transformation, reducing discontinuities at its edges. This can be done by multiplying the signal with a suitable window function.

Problem 4-2

a) Use the following Matlab code to apply rectangular, triangular and Hamming window. Observe the level of magnitude and the distribution of side-lobs carefully to find any improvement of leakage reduction. Put the screen display below and make your comments about different level of leakages. Make sure the shape of the signal in the time domain is drawn according to the shape of window function.

%Windowing for reducing leakage %sampling rate 1024Hz

N=1024; dt=1/N; t=0:dt:1-dt;

%****define the signal** x=sin(2*pi*72 .*t )+... sin(2*pi*196.5 .*t )+... sin(2*pi*408 .*t );

%*********************** % Apply rectangular window subplot(3,2,1); plot(t,x);

axis([0 1 -5 5]);

xlabel('nT (seconds)'); ylabel('x[n] rectang win'); X=fft(x); df=1;

f=0:df:N-1; subplot(3,2,2);

plot(f,20*log10(abs(X))); axis([0 N/2 -30 60]); xlabel('k (Hz)'); ylabel('|X[k]| (db)');

% Apply triangular window for n=1:N

w(n)=1-abs(2*n-N+1)/N; end;

x1=x .*w; subplot(3,2,3); plot(t,x1);

axis([0 1 -5 5]);

xlabel('nT (seconds)');

ylabel('x[n] triang win'); X=fft(x1); df=1;

f=0:df:N-1; subplot(3,2,4);

plot(f,20*log10(abs(X))); axis([0 N/2 -30 60]); xlabel('k (Hz)'); ylabel('|X[k]| (db)');

% Apply Hamming window for n=1:N

w(n)=0.54+0.46*cos((2*n-N+1)*pi/N); end;

x1=x .*w; subplot(3,2,5); plot(t,x1);

axis([0 1 -5 5]);

xlabel('nT (seconds)');

ylabel('x[n] Hamming win'); X=fft(x1); df=1;

f=0:df:N-1; subplot(3,2,6);

plot(f,20*log10(abs(X))); axis([0 N/2 -30 60]); xlabel('k (Hz)'); ylabel('|X[k]| (db)');

b)Use your code with Hamming window to analyse the signal:

252.6n2175.5n23.2n2365.3nx[n]sincossincos

1024102410241024

Plot the graphic result with appropriate title and labels.

5x[n] Hamming win0-500.10.20.30.40.50.6nT (seconds)0.70.80.916040|X[k]| (db)200-20050100150200250300k (Hz)350400450500