USA: +1(669)289-8683 
Sign in
Write a
MATLAB program for problem 16 at the
end
of Chapter 4. Answer all questions given in the
problem.
Your program must contain comments stating purpose of the program and explanations of
different parts of the program. Write a well organized and commented program.
At the start of the program, the program MUST use disp statements to state the program
purpose and how it can be terminated. Then it must get input from the user. If the
inputted Rs is not positive, then terminate the program.
Your program could use a while loop as follows.
…
Rs = input(‘Enter a resistor value: ‘)
;
while Rs > 0
…Rs = input(‘Enter a resistor value: ‘)end..
Provide a program listing, and a print out of problem (program) results.
Given is the circuit shown in Fig. P4.16.
(a) In terms of V, Rs and R, obtain an expression for the power P delivered to R.
(b) Write a script that receives as input the values of V and Rs, and makes a vector R with
element values that range from 0.1*Rsto 10*Rs in increments of 0.1*Rs. Use a for loop to
create a vector of power delivered to each element value of R. Then, plot the power P versus R.
Do this for V = 10 volts, Rs = 10 Ohms and Rs = 100 Ohms. For each case, from inspection of the
plots, give a value of R such that V delivers the maximum power to R. Provide a script listing
and plots, including axes labels and titles.
(c) Vectorize your expression for P, and repeat part(b) without using a for loop.
(d) In general, for a given value of Rs, what value of R enables the resistor R to receive the most
power from the voltage source V?
Write a
MATLAB program for problem 16 at the end of Chapter 4. Answer all questions given in the
problem.
Your program must contain comments stating purpose of the program and explanations of
different parts of the program. Write a well organized and commented program.
At the start of the program, the program MUST use disp statements to state the program
purpose and how it can be terminated. Then it must get input from the user. If the
inputted Rs is not positive, then terminate the program.
Your program could use a while loop as follows.
.
.
.
Rs = input(‘Enter a resistor value: ‘) ;
while Rs > 0
.
.
.
Rs = input(‘Enter a resistor value: ‘)
end
.
.
Provide a program listing, and a print out of problem (program) results.
ChapterFOUR
Roland Priemer, Copyright 2012, all rights reserved
___________________________________________________________________
Program Flow Control
A method to solve a problem may include alternative steps. Which next step to take may
depend on results from previous steps. A script that implements a solution method must be
able to follow alternative steps based on answers to questions like: is a variable value different
from an allowed set of values, is the result of a calculation zero, do two variables have the same
value, which of many variables has a particular value, and other questions. The answers to
these kinds of questions can influence what to do next in a script.
After you have completed this chapter, you will know how to:
• test conditions and execute alternative script segments
• repeatedly execute a script segment until a condition occurs
• execute one of many script segments depending on input data or intermediate results.
4.1 Relational Operators
Relational operators are used to compare variable values and produce a result that is true or
false. For example,
>> a = 2; b = 3; % assign values to two variables Table 4.1 List of Relational Operators
Relational Operator Description <= less than or equal to
> greater than
>= greater than or equal to
Table 4.1 gives the relational operators that can be used in MATLAB. Notice that the left_expression relational_operator right_expression The expressions on each side of a relational operator can be any MATLAB expressions that Example 4.1 ______________________________________________________ >> a = 2; b = 3; c = -1; d = 4; % define four arithmetic variables Suppose the variables a and b are real scalars. Then, one of the relations a < b
a == b
a > b
must be true. However, if a and b are replaced by A and B, which are matrices with the same Table 4.2 Functions That Return a Logical Value
Function Example, A and B are matrices, c is a logical scalar and C is a logical matrix C = any(A), if A is a matrix, then the function, any, works on each column of A, all c = all(A), are all elements of A nonzero, if A is a vector
Relational expressions evaluate to a logical value. Relational expressions can be 4.2 Logical Operators Logical operators work on logical variables and produce a logical result. For example, the MATLAB uses the symbols &, |, and ~ for the logical operators and, or and not, Table 4.3 Definition of Logical Operators
AND OR XOR NOT
Logical operators can work with numbers. If a number is 0, then it is used as if it is Example 4.2 ______________________________________________________ >> a = 0; b = 1; c = -3; % assign scalar values MATLAB also has built-in logical functions for the logical operators. i.e., and(a,b) = a&b, Logical operators can work with matrices when both operands have the same logical matrix, where a nonzero element of A results in a logic 1 element of ~(~A). The built-in Another useful built-in function is the function find to locate the elements of a matrix Example 4.3 ______________________________________________________ >> A = [1 -5 0; -5 0 4; 0 0 2] % assign a matrix Use parentheses to control the sequence in which logical operations are evaluated. Example 4.4 ______________________________________________________ >> a = 0; b = 1; c = 1; d = 0; % assign logical values 4.3 If – Elseif – Else – End With the if-elseif-else structure, a block of statements is executed only if a condition is met. if expression_1 block #1 of statements elseif expression_2 block #2 of statements insert as many additional elseif and block of statements as needed else block #N of statements end where each expression must evaluate to a logical value. When MATLAB encounters the 1) evaluate expression_1 It is useful to indent each block of statements to easily see that each if key word is matched by if expression block of statements Here, if the expression is true, then the block of statements is executed, and if it is false, then >> x = -6.8; Example 4.5 _________________________________________________________ Suppose we are given a data point x , where its value can be anywhere in the range ( ) /w b a N= − . If we know which bin contains x , then we can say that x is known to be ( ) (2 1) , 1, , 2
w
c k a k k N= + − = (4.1)
Suppose that all we know about x is that it is contained within a particular bin, say bin k . Instead of using bins that all have the same width, a bin center vector c with N There is another purpose of this script. It illustrates writing a script that is reasonably To find a bin number, the range [ , ]a b is normalized with
x a
y N −
=
−
Therefore, if x a= , then 0y = , and if x b= , then y N= , and the range of y is [0, ]N . For a % Program to find the bin number Program 4.1. Program to find the bin number of x in the range [a,b]. Notice that all statements within matched if, else, elseif, and end key words are indented for 4.4 For Loop With the for loop structure, you can repeat executing a block of statements a predetermined for loop_index = first_index : index_increment : last_index block of statements end Every for key word must be matched by an end key word. The block of statements can include For convenience, denote the loop_index variable by k, first_index by a, index_increment 1) check incr, and if incr is positive, check that a <= b, and if a <= b, then set k = a
2) execute the block of statements
3) go back to the for statement and set k = k + incr, which gives the new loop_index value
4) if k <= b, then execute the block of statements again and go back to step (3)
5) if k > b, then skip the block of statements, and continue execution after the end statement. If a>b, then the block of statements is not executed, and execution continues after the end If in step (1), incr is negative, then check that a >= b, and if a >= b, then set k = a, and do 2) execute the block of statements where v is a vector and k is first set to v(1) and the block of statements is executed . Then, If you do not include an index_increment, as in for k = a:b then a default index_increment value of 1 is used, which requires that a<=b.
Example 4.6 _______________________________________________________
Let us take a look at the details of multiplying two matrices with for loops. Let A be an
(NxK) matrix, and let B be a (KxM) matrix. Consider C = A*B with dimension (NxM). Prog. 4.2
finds C with for loops and the built-in MATLAB matrix multiply operation.
% Program to multiply two matrices (C=A*B) using for loops
clear all; clc
N=3; K=4; M=2; % specify matrix dimensions
A=[1 0 2 -1;2 1 -1 0; 1 -1 1 2]; % specify two matrices
B=[3 0; -1 5; 1 3; 2 1];
C = zeros(N,M); % preallocating space for C, which decreases execution time
t = tic; % tic, a built-in function, starts a stopwatch timer
for n = 1:N % index for rows of C
for m = 1:M % index for columns of C
sum = 0; % initializing the sum of products
for k = 1:K % get dot product of a row of A and a column of B
sum = sum + A(n,k)*B(k,m); % accumulate products
end
C(n,m) = sum; % save dot product of a row of A and a column of B in C
end
end
loops_compute_time = toc(t) % toc gets elapsed computing time since tic
C % display result
t = tic % start stopwatch timer
D=A*B; % check
MATLAB_compute_time = toc(t) % get elapsed computing time since tic
Program 4.2. Use of for loops to multiply two matrices. Notice that all statements within matched for and end key words are indented for for k = 1:K, sum = sum + A(n,k)*B(k,m); end This for loop could also have been replaced to find C(n,m) with
C(n,m) = A(n,:)*B(:,m); % dot product of row n in A and column m in B The script output is shown below. It was rearranged to take less space. loops_compute_time = 0.0024 secs This example demonstrates a programming utility of MATLAB, i.e., it works with the The built-in functions tic and toc were used in Example 4.6 to measure the execution If in a script the size of a matrix (or vector) is increased with each pass through a loop, it Example 4.7 ______________________________________________________ A signal is ( ) sin( )x t A tω= , and we want to see how the values of ( )x t are distributed ( ) (2 1) , 1, 2, , c k A k k N = − + − = (4.2)
Let 2ω π= rad/sec, and therefore one cycle takes 0 1T = sec. The signal ( )x t will be sampled at the % Program to find the histogram of a sine wave Program 4.3. Program to find the histogram of a sine wave. -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
0
500
1000
1500
2000
2500
3000
3500
4000
4500
sine wave value
va e s t
rib
ut n
sine wave histogram
Figure 4.1. Histogram of a sinusoidal signal. The histogram is shown in Fig. 4.1. We see that a sine wave spends almost equal amounts of This example introduced the built-in function bar, which plots a bar chart. Like the built- In the previous examples, the for loops always executed to completion. Sometimes it is Example 4.8 ______________________________________________________ Evaluate the sin( )x function using its Taylor series, which is given by
2 1
0 ( 1) k k
k x ∞ +
= −
In a script, the upper summation limit cannot be infinity. We must either set the upper % Program to evaluate the sin(x) function using its Taylor series % the size of y_approx is increased with each pass through this loop
Program 4.4. Evaluate the sin function using its Taylor series.
0 1 2 3 4 5 6 7 -1
– 0. 5
0 0.5
1 angle – radians
sin function with Taylor series
Figure 4.2. sin function produced with its Taylor series. The outer for loop in Prog. 4.4 iterates over the values in the vector theta. The inner for loop in max_k_value = 22
Some built-in functions concerned with controlling loop activity are given in Table 4.4. Table 4.4 MATLAB Functions to Exit For Loops, While Loops and User-Defined Functions
Function Brief Description skipped and execution continues with the next iteration of the loop in which break When break is encountered in a loop, it terminates execution of the loop. In return When return is encountered in a user-defined function, it causes an exit from
While loops will be discussed in the next section. 4.4.1 Probability To consider some interesting problems, a brief background in probability will now be given. Suppose we do an experiment, and the outcome is not predictable. Assume that an For example, consider tossing a coin. Each coin toss is an experiment, and there are two lim HH N p = (4.3)
and for tails, the probability Tp is given by
lim TT N p = (4.4)
Since H TN n n= + , we have a fundamental property of the probabilities of all possible 1lim limH T H TN Np p + = + =
For a fair coin we expect that 0.5Hp = and 0.5Tp = . If there are K possible outcomes of an experiment, and , 1, 2, ,kp k K= is the
1
1 k p =∑ (4.5) Practically, we cannot do an experiment an infinite number of times to find a probability We can use the MATLAB built in function rand to experimentally find probabilities x = rand; % 0 < x < 1
assigns to x a number in the range (0,1). Each use of rand returns a number that does not
depend on numbers obtained previously.
A pseudorandom number generator is an algorithm for generating a sequence of
numbers that approximates the properties of random numbers. The sequence is not truly
random. The function rand is a pseudorandom number generator. A careful mathematical
analysis is required to have any confidence that a pseudorandom number generator produces
numbers that are sufficiently random to meet the needs of a particular application.
Nevertheless, it is used in simulations to study the behavior of systems in which random events
can occur.
http://en.wikipedia.org/wiki/Algorithm� http://en.wikipedia.org/wiki/Random� The algorithm for the function rand is initialized the same way each time a MATLAB Example 4.9 _____________________________________________________ Let us simulate the coin toss experiment to find Hp and Tp . This will be done by clear all; clc; Program 4.5. Using a random number generator to simulate tossing a coin. The program output is shown below. probability of heads: pH = 0.4995 If we increase the number of tosses, then it is likely that pH and pT will each be closer to 0.5. The mean (average) value of the numbers produced by rand can be found rng(‘default’) % initialize random number generator x_mean = 0.5000 With the function rand we can obtain uniformly distributed random numbers over any Example 4.10 _____________________________________________________ To obtain random numbers x with a uniform distribution over the range: -5 < x < +5, use
>> N = 1e7; % specify the number of random numbers Let us see the first 100 points in x with >> y = x(1:100); % get first 100 points from x More about plotting will be given in Chapter 9. Fig. 4.3 shows a plot of some points from x. It 0 10 20 30 40 50 60 70 80 90 100 0 index
ra om ig l
Figure 4.3. A random signal. An important communications problem occurs due to noise. This is because a received ( ) ( ) ( )x t s t n t= + where ( )n t is additive noise. For example, the signal ( )x t could be the playback from an old ( )s t ( )x t
Process that operates on estimate of s(t)
Figure 4.4. Pictorial of a desired operation. To design an estimation process, it will be useful to do simulation studies with a computer to MATLAB also includes the two built-in functions randi and randn. The function randi matrix of random integers, where each element in A is in the range 1 to imax, while x = Before describing the random numbers produced by randn, let us study the random 2σ , of a set of random numbers is a measure of their volatility. The variance is computed with 2 2
1
1 ( ( ) ) n N =
where x denotes the mean of x , and σ is called the standard deviation. More specifically, The built-in function hist returns a vector, the elements of which are counts of the Example 4.11 _____________________________________________________ Prog. 4.6 finds the mean, variance and a histogram of the numbers obtained with % Histogram of numbers produced by the psuedo random number generator randn % this is the probability that a value of x will be in the range –sigma to Program 4.6. Script to find the mean, variance and a histogram of a set of numbers. The program output, which is given below, shows that the built-in function randn produces x_mean = -5.9235e-004 The function randn is designed to be a normally distributed pseudorandom number -4 -3 -2 -1 0 1 2 3 4
0
2 6
8
10 5 random number value
bi co t Distribution
Figure 4.5. Distribution of random numbers obtained with the function randn. If each bin count is divided by the total of number of samples N used to get the 0.125 0.125x− ≤ < is 0 0 / 10 5 / 1 7 0.1p b N e e= = . The function randn produces numbers in the range x−∞ < <+ ∞ , where a large value of abs( )x is very unlikely, and like (4.5) we have
that
1k p =−∞
=∑ (4.7)
______________________________________________________________ 4.4.2 Median Filtering In signal processing, it is often desirable to perform some kind of noise reduction on a signal The median filter concept is actually simple. Let us see how it works with a n=1, w=[5 5 6] → sort → u=[5 5 6], y(1)=u(2)=5. There is no x(0); x(1) was repeated. The outliers in x have been removed by median filtering. The signal y is less erratic (noisy) than The window w is said to slide across the signal x. For this demonstration K=3 was used. http://en.wikipedia.org/wiki/Signal_processing� http://en.wikipedia.org/wiki/Noise_reduction� http://en.wikipedia.org/wiki/Digital_filter� http://en.wikipedia.org/wiki/Digital_filter� http://en.wikipedia.org/wiki/Signal_noise� http://en.wikipedia.org/wiki/Image_processing� Example 4.12 _________________________________________________________ For an 88 key piano, let 1, 2, , 88n = be the key number from the lowest to the 1/12 ( 49)440 (2 ) nf −= Hz (4.8) For 40n = the note is middle C, and for 49n = , the note played is A440. Prog. 4.7 consists of two parts. In the first part, a user can enter an integer n , and then % Program to sound an 88-key piano Program 4.7 (first part). Sound a piano key. If the note produced in the first part of Prog. 4.7 comes from an old recording, then it noise and spike noise, which occurs randomly in the sound, and then it applies median filtering % example of median filtering % limit sound to +/- 1, required by audio device Program 4.7(second part). Median filter piano key sound. Fig. 4.6 shows note middle C with noise and spike noise, and the sound produced by this 1.135 1.14 1.145 1.15 1.155 1.16 1.165 -2
-1 3 seconds
am itu Figure 4.6. Note middle C with noise and spike noise. 1.135 1.14 1.145 1.15 1.155 1.16 1.165 3 seconds Figure 4.7. Median filtered note middle C. The signal shown in Fig. 4.7 still has some noise in it. However, additional digital signal 4.5 While Loop With the while loop structure you can repeat executing a block of statements until an while expression block of statements end Every while key word must be matched by an end key word. The block of statements can When MATLAB encounters the while key word it does the following 1) evaluate the expression A while loop can be made to work like a for loop. For example, n = 1; N = 10; end Within the block of statements, some condition must be tested that will either leave err Example 4.13 _________________________________________________________ The inner for loop in Prog. 4.4 of Example 4.8 could be replaced with k = 1; negligible = 0; Consider the problem of finding the value of x that minimizes a function ( )f x , assuming that We can see three minima in Fig. 4.8. At optx x= , ( )f x has its global minimum (also called the f(x)
xx opt x nx bx a
Figure 4.8. A function with a minimum at x = xopt.
Suppose nx is a guess of optx . In Fig. 4.8, nx is to the right of optx , where the slope of
( )f x at nx x= is positive, as indicated by the tangential line. Therefore, if the slope of ( )f x ‘1 ( )n n nx x f xµ+ = − (4.9) where µ , called the step size, is a small positive number that we must pick and ‘( )nf x is the To find a value of x that maximizes a function ( )f x , change (4.9) to ‘1 ( )n n nx x f xµ+ = + (4.10) The recursion, (4.9), of the method of steepest descent is operated for 0, 1, 2,n = , trial and error, and depends on the behavior of ( )f x as well as computation precision. If the the vicinity of a local minimum, such as ax or bx in Fig. 4.8, in which case (4.9) may never computing time) to converge. The recursion (4.9) (or (4.10)) is operated until ε+ − <1n nx x , for
some desired precision ε .
We do not have to keep the sequence nx , 0, 1, 2,n = . Equation (4.9) could be an
update of oldx to obtain newx , as in
1) set µ to a small positive number and set oldx to an initial guess of optx
2) compute ‘ ( )new old oldx x f xµ= − Example 4.14 _________________________________________________________ Find the value of t where 2te t− = . This is an example of a problem where the method minimizing 2 2( ) ( )tf t e t−= − , where ( ) 0f t ≥ for all t , and the minimum value is ( ) 0optf t = .
The square of 2( )te t− − is used to ensure that ( )f t is nonnegative. With
% Application of the method of steepest descent to solve the nonlinear end Program 4.8. Application of the method of steepest descent. f_of_t_initial_guess = 3.995764008937280e-001 The method of steepest descent has some drawbacks. It is slow to converge, which can
( / 2) ( / 2)
‘ ( ) n nn f x + − −
(4.11)
where h is some small positive number, for example, h µ= . Then, the method requires more The MATLAB built-in function fminsearch, which is a function function, is a derivative [x_opt,f_opt] = fminsearch(fun,x_initial_guess) where fun is a function handle. For example, fun could be @myfun, where myfun is the name Let us apply fminsearch to the problem of Example 4.14. An anonymous function is >> format long e
While fminsearch is easy to use and can give good results, it can also require many 4.7 Numerical Integration Given is some function ( )x t that is finite for t−∞ < < + ∞ . To find ( )y t given by
( ) ( ) ( ) ( ) , a y t x d x d x d t aτ τ τ τ τ τ = = + >∫ ∫ ∫ (4.12) is a common problem. We can write (4.12) as ( ) ( ) ( ), a y t x d y a t aτ τ= + ≥∫ (4.13) where ( )y a is called the integration constant. If t a= , then (4.13) becomes ( ) ( )y t a y a= = . More specifically and for example, recall Fig. 2.2, in Chapter 2, where
( )
( ) i t =
If we integrate this equation, we get ( ) ( ) ( ) ( ) ( ) ( ) , a a
q t i d i d i d i d q a t aτ τ τ τ τ τ τ τ = = + = + ≥∫ ∫ ∫ ∫ (4.14) The second term ( )q a is the amount of charge that has passed through the wire cross-section Depending on ( )x t in (4.13), it may be difficult, if not impossible, to find a function for
( )y t , and we must resort to finding ( )y t for a set of discrete time points with some numerical integration method. Sometimes we do not know ( )x t for all t . Instead, we only know ( )x t for Without loss of generality, let 0a = , and set the upper integration limit to 0t T= .
0
0 ( ) A x t dt= ∫ (4.15) where A is the area under the curve ( )x t from 0t = until 0t T= . 4.7.1 Euler’s Method Let us segment the time range 0[ 0, ]T into N sub-intervals, where the width of each sub- 1 1
0 0
( ) ( ) N N
k k − − = =
=∑ ∑ (4.16) where (( 1) )x N T− is the value of ( )x t at the beginning of the last sub-interval. The ( )x t tT 2T NT (0)x (2 )x T 0( )x T NT=
(4 )x T Figure 4.9. Segmentation of ( )x t over 0[0, ]T NT= . In view of (4.16), let 1 0 n k − = ∑ ( 1) ( ) ( ), 0, 1, , 1, (0) 0S n S n x n T n N S+ = + = − = (4.17) This is a recursion in ( )S n , where (1) (0) (0 ), set S(0) = 0 (3) (2) (2 )
( ) ( 1) (( 1) )
S S x T
x T
S S x T S N S N x N T
= +
= + = − + − Each time a new value ( )x nT becomes available, ( )S n is updated to obtain ( 1)S n + . The Example 4.15 _____________________________________________________ Let ( )x t be the function defined by
1, 0
( ) sin( ) t t π
π = = (4.18)
This function is called the sinc function, i.e., ( ) sin c( )x t t= , and it is a built-in MATLAB function. >> t = linspace(-10,10,1000); % specify time points we get the plot shown in Fig. 4.10. The sinc function oscillates indefinitely in both directions -10 -8 -6 -4 -2 0 2 4 6 8 10 0 t – secs
si (t Figure 4.10. The sinc function plotted from t=-10 to t=10. An accurate value for the integral of the sinc function from 0t = until 10t = is: % Program to find the integral of a function using Euler’s method Program 4.9. Program to find the integral of sinc( )t using Euler’s method. N Area These results show that N must be large to achieve a small numerical integration error. ______________________________________________________________ 4.7.2 Trapezoidal Rule The trapezoidal rule is a slight modification of Euler’s method, and gives better results. Instead 1 1 ( ) (( 1) ) 2 2 ( 1) ( ) ( ) (( 1) ), 0, 1, , 1, (0) 0
N N x kT x k T T S n S n x n T x n T n N S
− − + + + = + + + = − =
∑ ∑
(4.19)
where ( / 2) ( )A T S N . Example 4.16 _____________________________________________________ To apply the trapezoidal rule to the integration problem of Example 4.15, the inner for xnT = sinc(0); % this is x(nT) for n=0 and it will be updated in the loop and the area is computed with Area(m) = (T/2)*S; The results are: N Area Here, the result for 50N = are better than the result for 550N = using Euler’s method. 4.7.3 Built-In Integration Functions MATLAB has several sophisticated built-in functions for evaluating a definite integral. One of area = quad(@my_function,a,b, tol) where the name of the m-file that defines the function to be integrated from a to b must be The function quad was used to obtain the accurate result from integrating the sinc When MATLAB executes the quad function algorithm, it will want to pass to the user- Example 4.17 _____________________________________________________ Let us find the amount of charge that has passed through a diode over one cycle of a function current = diode(t) Here, the function input variable is t , which is used as the input time argument for the % Program to plot charge movement through a diode over one cycle of frequency = 60; % specify frequency Program 4.10. Script to find the amount of charge that has passed through a diode. 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
time secs
Co m Charge movement through a diode
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.8
– 0.6
– 0.4
– 0.2
0 Figure 4.11. Charge movement through a diode over one voltage cycle. Fig. 4.11 shows how charge moves through a diode over one cycle of the voltage across There are several other built-in functions for numerical integration. These functions may be more efficient than quad, while also giving more accurate results, especially for 4.8 Switch – Case – Otherwise The switch-case-otherwise structure provides a way to select which block of statements in a switch expression block #1 of statements case expression_2 block #2 of statements insert as many additional case expression and block of statements as needed otherwise block #N of statements end where each switch key word must be matched by an end key word. The otherwise and block 1) evaluate the switch expression and
It is useful to indent each block of statements to easily see each case block of statements and Example 4.18 _____________________________________________________ You are approaching an intersection with a traffic light. Different actions are required % Suggested action when approaching an intersection with traffic lights Program 4.11. Program to suggest action when approaching an intersection with traffic lights.
______________________________________________________________ 4.9 Conclusion Testing a condition, executing alternative parts of a script and repeatedly executing a part of a • to use relational operators for formulating condition tests • the if-elesif-else structure is used to execute alternative script parts on a condition
Table 4.5 gives the MATLAB functions that were introduced in this chapter. Use Table 4.5 Built-In MATLAB Functions Introduced in This Chapter
Function Brief Explanation and contain the same values, and logical 0, meaning false, L = and(A,B) an element of the output array is set to 1 if both input arrays L = or(A,B) an element of the output array is set to 1 if either input array L = xor(A,B) an element of the output array is set to 1 if either input array L = not(A) an element of the output array is set to 1 if A contains a zero && short logical and the array A it equals zero and -1 if it is less than zero; for the non-zero tic start a stopwatch timer T = toc saves the elapsed time in T (seconds) since the last execution of pause(t) pauses for t seconds before continuing, where t can also be a bar(X,Y,width) draws the columns of the M × N matrix Y as M groups of N vertical continue pass control to the next iteration of for or while loop keyboard; also terminates the keyboard mode control to the user’s keyboard. The special status is indicated by a R = rand(N,M) returns an N× M matrix containing pseudorandom numbers I = randi(N,M,imax) returns an N× M matrix containing pseudorandom integer values R = randn(N,M) returns an N× M matrix containing pseudorandom values from Y = mean(X) Y is a row vector containing the mean of each column of X sound device; X can be a two column matrix for stereo. Clips soundsc(X,fs) same as sound, but autoscales file.wav, with a sampling rate fs Hz; stereo data must be specified Y = sinc(X) returns in the matrix Y elements found with sin(pi*x)/(pi*x) of the quad(func,a,b) returns an approximation of the integral from a to b of the plotyy(x1,y1,x2,y2,plot_func) plots y1 versus x1 with y-axis labeling on the left and plots y2 quadl same as quad, but uses a more accurate numerical integration quad2d(func,a,b,c,d) returns an approximation of the double integral of a function ≤ ≤c y d . See MATLAB help for more options. The essentials of the MATLAB programming language have been presented in the first In the next chapter the concept of a logical variable will be extended to understand its References
[1] Nelder, J.A. and Mead, R. (1965), “A simplex method for function minimization”, Comput. J.,
[2] Li, C., Priemer, R. and Cheng, K.H. (2004), “Optimization by Random Search with Jumps”, Int. Problems In all of the scripts that you are required to write for the following problems, include comments that Section 4.1 1) Use a = -2, b = 3 and c = 5, and find the result of each of the following relational expressions. (a)
6) For A = [0 1; 1 0], C = [0 0; 1 0] and C = A & B, give all possible B. 7) For a = -2, b = 4 and c = 2, give the results of: (a) d = (a ~= c -2) & a^2 > b, 8) For A = [-1 0.5 5; 2 -4 0; 0 0 1], (a) give y = find(A <= -1 | A >= 1) , and explain the meaning of Section 4.3 9) (a) If the inputs to Prog. 4.1 are: a = 0, b = 5, N = 9 and x = 5.6, then what is the script output? 10) (a) Convert Prog. 4.1 into a function, call it bin_number. The function input arguments are: a, 11) Write a script that implements the flowchart given in Fig. P4.11. Use an if-else-end control Input a resistor value R
R > 0 ? Terminating Program Input a current I
Calculate power Display End program
Figure P4.11. Flowchart to calculate power delivered to a resistor.
12) Write a script that inputs an integer N in the range 0 to 9. Your script must check that the Section 4.4 13) Describe what is wrong with each of the following for statements. (a) for k = 0:-1:10, (b) for n = 10:2:3, (c) for 2k = 0:5, (d) for Number = -1:-5. In each case, 14) (a) Write a function script, call it inner_prod, that uses a for loop to find the inner product 0ω . For each K, what is A_dot_B? What do you think A_dot_B will be for any integer K ≠ 1?
(c) Repeat part (b) for B = sin(K 0ω t). 15) Write a script that uses for loops to evaluate the power series given by (2.13) of an
exponential function 2( ) tx t e−= , for t = 0:0.1:2. Break out of the for loop that does the power 16) Given is the circuit shown in Fig. P4.16. V R
R s
Figure P4.16. Battery with internal resistance Rs connected to a load R.
17) Write a script that starts by obtaining a vector x of N = 1e6 normally distributed random 18) Suppose that an experiment consists of tossing two die. Each die will show an integer from 1 to 19) Let us study the performance of the algorithm given by
1
( ) ( ( ) 2 (( 1) ) (( 2) ) y k T x k T x k T x k T= + − + −
for 0, 1,k = , where ( )x k T comes from sampling a sound given by and 2 2 (440)ω π= . The sampling rate is 8000sf = samples/second, and 1 / sT f= . The % This is an example of digital signal processing % since there is no x(-T) use Section 4.5 20) Do Prob. 4.11 using a while loop. The script should allow for repeated execution of the script clear all; clc; Explain why the abs function must be used. Complete the script, and provide results as
Section 4.6 22) Consider again Prob. 4.16. Write a script that applies the method of steepest ascent, (4.10), to 23) Given a number x, find y = x1/2. Of course, y = sqrt(x). However, an algorithm must be 24) Repeat Prob. 4.23, but replace the inner while loop with an invocation of the MATLAB built-in Section 4.7 25) Given is the function x(t) = e-at – e-bt. Let a = 1 and b = 2. 26) Repeat Prob. 4.25, but use the trapezoidal rule. 27) For the x(t) given in Prob. 4.25, use the built-in MATLAB function quad to find the area A A found in Prob. 4.25, part(a). Use format long e for all script output. Provide a listing of Section 4.8 28) Write a function script, call it linear_eqs_solver, where the inputs are a square matrix A, a 29) Repeat Prob. 4.18, where within the for loop a switch-case-end structure is used to count the
>> c = a < b % compare a and b using the relational operator <, meaning less than
c = 1
Since the value of a is less than the value of b, the comparison, a
< less than
== equal to
~= not equal to
relational operator, ==, which tests equality, is not an assignment operator. The syntax of a
relational expression is
evaluate to scalars, vectors or matrices having the same dimension. A relational expression
produces a logical result, which can be assigned to a variable.
>> e = sqrt(a^2+b^2) <= sqrt(c^2+d^2) % arithmetic expressions on each side of the relational operator
e = 1
>> f = (a > b) +2* (c < d) % logical values (and variables) can be used in arithmetic expressions
f = 2
>> % f, an arithmetic variable, can be 0, 1, 2 or 3
>> w = 2*pi; t = 0:0.1:1; % assign a frequency and a time range
>> X = sin(w*t) <= cos(w*t) % compare two vectors element-by-element, X is a logical vector
X = 1 1 0 0 0 0 0 1 1 1 1
>> % matrices that have the same dimensions can be used on each side of a relational operator
>> g = exp(c*t) .* sin(w*t); % define an arithmetic vector, an exponentially decaying sinusoid
>> Z = g > 0 % check times at which exponentially decaying sinusoid is positive
Z = 0 1 1 1 1 1 0 0 0 0 0
% Z is a logical vector, resulting from comparing each element of g to zero
>> g_pos = Z .* g % get positive values of g, a logical vector can be used in arithmetic expressions
g =
Columns 1 through 11
0 0.5319 0.7787 0.7046 0.3940 0.0000 0 0 0 0 0
_____________________________________________________________
dimensions, then the relation is applied to corresponding elements of A and B, and the result is
a logical matrix with the dimension of A, where some elements are logic 1, while others are
logic 0. Since A and B are matrices, then these relations do not produce a scalar logical value.
To test the equality of two matrices, use the built-in MATLAB function isequal, as in c =
isequal(A,B), where c is a scalar logical variable. Table 4.2 gives some built-in MATLAB functions
that operate on matrices and return a logical value.
isequal c = isequal(A,B), are all corresponding elements equal
islogical c = islogical(A), are all elements either 0 or 1
isinteger c = isinteger(A), are all elements integers
isfloat c = isfloat(A), are all elements floating point numbers
isempty c = isempty(A), does A have no elements
true C = true(M,N), C is an MxN logical matrix of all logic 1 values
false C = false(M,N), C is an MxN logical matrix of all logic 0 values
logical C = logical(A), convert numeric matrix A to logical matrix C
any c = any(A), is any element of A nonzero or logical 1, if A is a vector
and C is a vector. To obtain a scalar, use c = any(any(A))
C = all(A), if A is a matrix, then the function, all, works on each column of A,
and C is a vector. To obtain a scalar, use c = all(all(A))
compounded with the logical operators and, or, xor and not.
expression, (a < b) and (c == d), can only be true if both relations are true. It is false if either
relation or both relations are false. In MATLAB the word, "and", is replaced by the symbol &,
and a MATLAB assignment statement becomes e = (a
respectively, and a built-in function xor for the logical operator xor (exclusive or). Table 4.3
defines the logical operators, where a, b, and c are scalar logical variables.
a b c = a & b a b c = a | b a b c = xor(a,b) a ~a
0 0 0 0 0 0 0 0 0 0 1
0 1 0 0 1 1 0 1 1
1 0 0 1 0 1 1 0 1 1 0
1 1 1 1 1 1 1 1 0
logical 0, and if it is nonzero, then it is used as if it is logical 1.
>> a | b % both operands have logical values
ans = 1
>> b & c % c does not have a logical value, but it is used as if it is logic 1
ans = 1
>> ~c
ans = 0
>> xor(a,c) % produces a logic 1 only if one of the two inputs is logic 1
ans = 1
______________________________________________________________
or(a,b) = a|b, and not(a) = ~a. There is an option when using the & and | operators. Consider
c=expression_a & expression_b. If expression_a evaluates to logic 0 then there is no need to
evaluate expression_b, since c=0. To require that MATLAB evaluates expression_b only if
expression_a is logic 1, use && instead of &. Similarly, consider c=expression_a | expression_b.
If expression_a is logic 1, then there is no need to evaluate expression_b, since c=1. To require
that MATLAB evaluates expression_b only if expression_a is logic 0, use || instead of |. These
options can reduce execution time.
dimension. The logical operation is done element by element, and the result is a logical matrix.
If one operand is a scalar, then the elements of the other operand, a matrix, are used element
by element to produce a logical matrix. If A is a matrix of numbers, then ~A produces a logical
matrix where a nonzero element of A results in a logic 0 element in ~A. And, ~(~A) produces a
functions given in Table 4.2 can be applied to logical matrices.
that are nonzero. The function returns the indices numbered linearly down each column
starting at the left column.
A =
1 -5 0
-5 0 4
0 0 2
>> % the indices for the first column are: 1, 2 and 3, and for the second column
>> % the indices are: 4, 5 and 6, and so on
>> find(A) % find the linear indices of all elements in A that are nonzero
ans = 1 2 4 8 9
>> B = [0 -1.8 -4 5.2 0.6 0 1];
>> I = find(B) % find the indices of the nonzero elements of a vector
I = 2 3 4 5 7
>> find(B,3,’first’) % find the first 3 indices of the nonzero elements in a vector
ans = 2 3 4
>> find(B,4,’last’) % find the last 4 indices of the nonzero elements in a vector
ans = 3 4 5 7
% logical operations and relational operations can be combined in compound expressions
>> find((B > -1) & (B < 1)) % find the indices of the elements in the range: -1 < element < 1
ans = 1 5 6
% Here, (B > -1) produces a logical vector, which is anded element by element with the logic vector
% produced by (B<1) to produce the logical vector that is the input to the find function
______________________________________________________________
MATLAB evaluates an expression within the innermost parentheses first.
>> a & b | c & ~d % MATLAB evaluates the logical and operations first
ans = 1
% this results from: (a&b) | (c & ~d)
>> a & (b | c) & ~d % MATLAB evaluates the term within parentheses first, and then left to right
ans = 0
% parentheses must be carefully used to achieve the desired overall logical operation
___________________________________________________________________
The syntax of the most general form of this structure is
if key word, it does the following
2) if expression_1 is true, then execute block #1 statements, and continue execution after the
end statement
3) if expression_1 is not true, then skip the block #1 statements to the first elseif key word and
evaluate expression_2
4) if expression_2 is true, then execute block #2 statements, and continue execution after the
end statement
5) if expression_2 is not true, then skip the block #2 statements to the next elseif key word, if
there is one, evaluate its expression, and so on
6) if all expressions are false, then MATLAB will skip to the else key word, if there is one, and
execute the block #N statements, after which execution continues after the end statement.
an end key word. The elseif through the else sections are optional. Within each block of
statements there can be nested if-elseif-else-end structures. You can nest these structures as
deeply as needed. The simplest if structure is
end
the block of statements is skipped and execution continues after the end statement. If the
block of statements can fit on one line, then the if statement could be, for example
>> % limit x to the range: -5 to +5
if (x<-5) | (x>5), x = 5*sign(x); end % the built-in function sign returns -1, 0 or +1
>> x
x = -5
from a to b, written as [ , ]a b . The brackets indicate that the end point values are included in
this range. If we write ( , )a b , then both end points are not included in this range. It is desired
to know which of N sub-intervals of the range contains x . Let us segment this range into N
segments (sub-intervals). Each segment is called a bin. The width of each bin is
within its bin range. The centers of the bins are given by
Then, given k only, the best estimate of the value of x is ( )c k , and the maximum estimation
error is / 2w .
elements could be given, where the bin ranges are: [ , ( (1) (2)) / 2]a c c+ ,
[( (1) (2)) / 2, ( (2) (3)) / 2]c c c c+ + , , [( ( 1) ( )) / 2, ]c N c N b− + , for a total of N bins over [a,b]
that have varying widths. Using bins that all have the same width, Prog. 4.1 finds 1, ,k N= ,
the bin number, given x .
robust. This means that input data or intermediate results within the script do not cause a
computer to become non-responsive or generate error messages. A script may be robust and
not give correct (or expected) results, because its programming is based on an incorrect
algorithm. If a script executes to completion and the results are incorrect, then the
programmer (possibly working with a user), must interpret the results, and make corrections.
b a
given x , ceil( )k y= , where ceil is a built-in function that returns the next higher integer.
clear all; clc
disp(‘Program to find the bin number of x in the range [a,b]’)
err = 0; % used for checking if an error occurred
a = input(‘enter lower range limit: ‘);
b = input(‘enter upper range limit: ‘);
if a < b % a robust program must check inputs
N = input('enter the number of bins, greater than 0: '); N=round(N);
if N > 0 % check input
x = input(‘enter x: ‘);
if x >= a && x <= b % check input
w = (b-a)/N; % get bin width
k = ceil((x-a)/w); % round to next higher integer
else
err=1; % x is out of range
end
else
err=2; % N is out of range
end
else
err=3; % lower limit not less than upper limit
end
if ~err % was an error detected
fprintf('out of N = %i bins, the bin number is: %i \n',N,k)
elseif err == 1, disp('x is out of range')
elseif err == 2, disp('N is out of range')
elseif err == 3, disp('lower limit not less than upper limit')
end
clarity.
___________________________________________________________________
number of times. Its syntax is
additional for loops. For loops can be nested. The block of statements can include if-elseif-
else-end structures and any other valid MATLAB statements. It is useful to indent the block of
statements to easily see the matched for and end key words. The loop_index must be a
MATLAB variable. The first_index value, index_increment value and last_index value can be
expressions that evaluate to a real scalar.
by incr, and last_index by b. When MATLAB encounters the for key word it does the following
statement.
the following
3) go back to the for statement and set k = k + incr, which gives the new loop_index value
4) if k >= b, then execute the block of statements again and go back to step (3)
5) if k
k=v(2) and the block of statements is executed, and this continues until all elements of v have
been used. With this structure you can assign to k any values, including complex values, for
example, v = [1-j*2 2* pi -1 exp(j*pi/4) u], where u = [1/3 j*3].
block of statements
end
clarity. Since the innermost loop has only a few statements, it could have been written as
C =
3 5
4 2
9 0
MATLAB_compute_time = 1.4848e-005 secs
matrix data type. Moreover, D=A*B is executed in much less computing time than it took to
obtain C, even though space for C was preallocated.
___________________________________________________________________
times of parts of Prog. 4.2. The time required to execute a part of a script depends on the
speed of your computer. There may be occasions when, for example, between different
program outputs, you want to control the time in between these program outputs. To do this,
use the built-in function pause(seconds), for example, pause(4.5), which causes a program to
stop execution for 4.5 secs. Use pause without an argument to stop a program and wait for the
user to strike any key before continuing.
is best to preallocate space for the largest size of the matrix before getting into the loop. This
will save much execution time.
over one cycle. The signal value varies over the range A− to A , written as [ , ]A A− . Let us
segment this range into N bins. The width of each bin is 2 /A N , and using (4.1) the bin
centers are given by
A
N
rate 44,100sf = samples/sec. Prog. 4.3 counts the number of times that the value of ( )x t is in
each of the N bins. The resulting distribution is called a histogram.
clear all; clc
% specify the data
A = 1; freq = 1; T0 = 1/freq; % specify amplitude and freuency
fs = 44100; T = 1/fs; % sampling rate and time increment
t = 0.0:T:T0-T; % time points over one cycle
x = A*sin(2*pi*freq*t); N_data = length(x); % get samples of x(t)
N = 40; w = 2*A/N; k = 1:N; % number of bins, bin width, and bin index
c = -A+(2*k-1)*w/2; % bin centers
bin = zeros(1,N); % preallocate space and initialize bin counts to zero
for n=1:N_data
k = ceil((x(n)+1)/w);
if k==0
k=1; % x=-1 belongs to bin 1
end
bin(k)=bin(k)+1; % increment the bin count
end
bar(c,bin,1) % get a barchart plot
grid on
% the label were entered using the plot editor
lu
di
io
time in segments from -0.8 to +0.8, while nearly 20% of the time in a cycle is spent within 5% of
the maximum and minimum values.
in function plot, the first input gives horizontal plotting information, the bar centers, the second
input gives the bar heights, and the third input controls spacing between bars in the plot, where
a value less than 1 causes a space between bars. More detail will be given in Chapter 9.
______________________________________________________________
useful to terminate a for loop before it has iterated through its entire index range.
sin( )
(2 1)!
x
k
=
+∑
limit to a large integer or use some way of stopping the computation after some predetermined
precision has been achieved. When x is incremented by 2π the function sin( )x repeats itself.
Let us evaluate the function for 0 2x π≤ ≤ . Program 4.4 evaluates the Taylor series and plots
the result. The resulting approximation of the sin function is shown in Fig. 4.2.
clear all; clc
N = 101; % N values for x will be used
theta = linspace(0,2*pi,N); % assign values to theta
y = sin(theta); % use the built-in function sin to compare with the series
K = 100; % K terms in the Taylor series will be used
precision = 1e-6; % this is a precision of 0.0001 percent
y_approx = zeros(1,N); k_value = zeros(1,N); % preallocate space
for n=1:N % loop for each value in theta
x = theta(n);
% no need to evaluate the Taylor series for x=0
if x == 0.0, y_approx(n) = 0; k_value(n) = 0; continue; end
sum_T = 0; % initialize Taylor series sum
factorial_term = 1; % initialize (2k+1)!
minus_one_term = 1; % initialize (-1)^k
x_power_term = x; % initialize x^(2k+1)
x_square = x^2; % used to get next x^(2k+1), computed only once
next_term = minus_one_term*x_power_term/factorial_term; % k=0 term
for k=1:K % loop to accumulate series terms
sum_T = sum_T + next_term; % accumulate terms
% prepare the next term in the Taylor series
factorial_term = factorial_term*(2*k)*(2*k+1); % update (2k+1)!
% uses only two multiplies to find each (2k+1)!
minus_one_term = – minus_one_term; % update (-1)^k
x_power_term = x_power_term*x_square; % update x^(2k+1)
% uses only one multiply to find each x^(2k+1)
next_term = minus_one_term*x_power_term/factorial_term;
% compare next term in Taylor series to sum and check precision
if abs(next_term) >= abs(precision*sum_T);
continue % continue this for loop
else
sum_T = sum_T + next_term; % include the next term
break % jump out of this for loop
end
end
% space for its largest size was preallocated before this loop
y_approx(n) = sum_T;
k_value(n) = k; % save index for which this for loop achieved precision
end
% finished evaluating the Taylor series for each element of theta
max_k_value = max(k_value) % get maximum number of for loop iterations
max_error = max(abs(y-y_approx)) % get maximum error
plot(theta,y_approx)
grid on
xlabel(‘angle – radians’)
title(‘sin function with Taylor series’)
Prog. 4.4 is set up to iterate K=100 times. It computes each term in the Taylor series. Since
K=100 terms may not be necessary to achieve a desired precision, each new term is checked to
see if it has become small enough to make a negligible contribution to the summation. If it has
not become small enough, the built-in function continue is used to continue the inner for loop
iteration. If it has become small enough, then the built-in function break is used to exit the
inner for loop in which break occurs. Some script parameters of interest are shown below.
While the inner for loop could have iterated up to K=100 times, the maximum number of
iterations to achieve the prescribed precision of 0.0001 % over a cycle of the sin function
turned out to be k=22, and the resulting maximum error when compared to sin( )x using the
MATLAB built-in function is 3.7e-6 percent. The use of tic and toc will show that the algorithm
that MATLAB uses to evaluate the sin function is much faster than a direct application of its
Taylor series.
max_error = 3.7065e-008
______________________________________________________________
continue When continue is encountered in a loop, remaining loop statements are
it appears. In nested loops, continue causes the next iteration of the loop
enclosing it.
nested loops, break cause an exit from the loop enclosing it only.
the function. Statements in the function after the return statement are not
executed. Also, terminates the keyboard mode, which will be discussed
Chapter 10.
outcome can be one particular event among a finite number of possible events. The probability
of an outcome gives an assessment of how likely it is that a particular outcome will occur
compared to all other possible outcomes.
possible outcomes or events, i.e., heads or tails. Suppose you toss a coin N times. Of the
N tosses, let Hn be the count of the outcome heads, and let Tn be the count of the outcome
tails. The probability Hp of heads is given by
n
N→ ∞
n
N→ ∞
outcomes, which in this case is:
n n
N N→∞ →∞
probability of each outcome, then
K
k
=
by counting. However, if we make N large enough, then we might get useful values for the
probabilities. The probability of an impossible outcome is zero, while the probability of a
certain outcome is one. For the probability p of any outcome we have 0 1p≤ ≤ . Notice that if
an outcome occurs only once in an infinite number of trials of an experiment, then its
probability is zero. This means that if the probability of an outcome is zero, we cannot be
certain that it will never occur. Similarly, if the probability of an outcome is one, then we
cannot be certain that it did not occur at least once.
associated with the outcomes of some experiments. In MATLAB, each rand function evaluation
returns a number in the range: 0 1→ (written as (0,1)), which we cannot predict, and where
all numbers in this range are equally likely, which means that the function rand is a uniformly
distributed random number generator. The MATLAB statement
session is started. Therefore, when you start a MATLAB session, rand will produce the same
random sequence of numbers that it produced the previous time that a MATLAB session was
started. During a MATLAB session you can reset the random number generator to the default
MATLAB start-up state with the statement: rng(‘default’). Use the MATLAB help facility to find
out more about the rand function.
counting the number of times each outcome occurs over a large number of coin tosses. Prog.
4.5 simulates the coin toss experiment.
rng(‘default’); % initialize random number generator to its default state
total_tosses = 1e6; % specify the total number of coin tosses
heads = 0; tails = 0; % initialize outcome counters
for n = 1:total_tosses
% do experiment
x = rand; % get a random number distributed from 0 to 1
if x < 0.5 % let this result be analogous to heads occurring
heads = heads + 1;
else % if not heads, then tails
tails = tails + 1;
end
end
disp('probability of heads: ');
pH = heads/total_tosses
disp('probability of tails: ');
pT = tails/total_tosses
probability of tails: pT = 0.5005
However, this does not mean that as the total number of tosses N is increased, then the
number Hn of heads or the number Tn of tails will ever become exactly / 2N . Since pH is
essentially the same as pT, the outcomes heads and tails are equally likely.
_______________________________________________________________
experimentally with
N = 1e7; % specify the number of random numbers to get
% obtain a vector of N random numbers uniformly distributed from 0 to 1
x = rand(1,N);
x_mean = sum(x)/N % find the mean value
% this can be obtained with the MATLAB function mean to get x_mean=mean(x)
range (a,b). The span of this range is (b-a). Since the range of x = rand is (0,1), the range and
span of y = (b-a)x are (0,(b-a)) and (b-a), respectively. Therefore, the range of y = (b-a)x + a is
(a,b).
>> x = rand(1,N); % get a vector of N uniformly distributed random numbers between 0 and 1
>> x = 10*x – 5; % multiply by (5-(-5))=10 and translate by -5
>> x_mean = mean(x) % find the mean value
x_mean = 2.9411e-006
>> % ideally, the mean value of x should be x_mean = 0.
>>
>> y = y – mean(y); % make the mean of y zero
>> plot(y,’- .’); grid on; xlabel(‘index’); ylabel(‘random signal’)
>> % within single quotes, the dash causes the plot function to connect the points with straight
>> % lines, and the period causes the plot function to put a dot at each point
looks like a noise signal with a zero mean value.
-5
5
nd
s
na
______________________________________________________________
signal ( )x t is a noise corrupted version of a transmitted signal ( )s t , such as
audio recording, and ( )s t is the audio signal that was originally recorded. Or, ( )x t could be a
received cell phone signal, and ( )s t is the signal that was transmitted from a distant location.
We are interested to design a process that performs as depicted in Fig. 4.4, where the process
output ˆ( )s t is an estimate of ( )s t .
its input x(t) to produce an
test the performance of various process designs and noise signals. This will require simulating
the noise signal with a model of noise, such as the noise signal shown in Fig. 4.3. Such
simulations are often done with MATLAB.
produces uniformly distributed random integers. For example, A = randi(imax,M,N) is an (MxN)
randi(imax) is a scalar integer in the range 1 to imax.
numbers it produces by finding the mean, variance and a histogram. The variance, denoted by
N
x n x
σ
= −∑ (4.6)
the variance (and standard deviation) of a set of random numbers is a measure of the degree to
which the random numbers in the set deviate away from their mean. MATLAB has a built-in
function, var, that computes the variance. A histogram of a set of random numbers shows how
the numbers are distributed.
number of times the elements of the vector of random numbers have values within the range
of each of the bins. We could have used the hist function in Example 4.7.
randn.
clear all; clc
N = 1e7; % specify the number of random numbers to get
x = randn(1,N); % obtain a vector of N random numbers
x_mean = mean(x) % find the mean value
x_var = var(x) % get the variance of the random numbers in x
bin_width = 0.25; % specify a bin width less than 1
N_pos_x_bins = 25; % specify number of bins for positive x(n) values
pos_x_bin_center = N_pos_x_bins*bin_width; % right most bin center
% locate equally spaced bin centers
bin_centers = -pos_x_bin_center: bin_width : pos_x_bin_center;
fprintf(‘using %i bins \n’,length(bin_centers))
x_hist = hist(x,bin_centers); % return a vector with the count in each bin
% prob_x = x_hist/N; % probability of x within each bin
bar(bin_centers,x_hist,1,’w’); % full bar width; and white color
axis([-4 4 0 1.05*max(x_hist)]); grid on
xlabel(‘random number value’); ylabel(‘bin count’)
title(‘Distribution’);
bins_to_1 = round(1/bin_width); % number of bins accounting for 0
% +sigma, where sigma is the standard deviation. With randn, sigma = 1
% additional plot editing was done in the Figure Window
random numbers having a zero mean and a variance of one. The standard deviation is one, and
within ± one standard deviation 67.92 % of the returned values occur. The histogram has a bell
shaped appearance. A normally (also called Gaussian) distributed random number sequence
has these properties.
x_var = 0.9993
using 51 bins
prob_minus_to_plus_sigma = 0.6792
generator, and produces numbers with zero mean and unity variance. With y a x b= + you
can produce a normally distributed random number sequence y having a prescribed mean and
variance, where 2 2y aσ = and y b= .
4
x 10
n
un
histogram of Fig. 4.5, then, like (4.3) and (4.4) , we get the probability that a number will occur
in a particular bin. Designate the bin counts for positive (negative) numbers by kb ( kb− ),
respectively. Then, for example, 0b is the bin count for 0.125 0.125x− ≤ < , and the probability
that an experiment outcome (invocation of randn) will result in a value of x within the range
k
+∞
(one-dimensional or 1-D signal) or an image (2-D signal). Median filtering is a digital signal
processing algorithm that is easy to apply and often used to reduce noise. For example, an old
audio recording, where there is crackling noise, can be improved by median filtering. Median
filtering of a speech signal is sometimes a pre-processing step before speech recognition
algorithms are applied. Median filtering is widely used in digital image processing.
demonstration. Suppose a signal s(t) is sampled, like in (3.33), to produce the number
sequence (discrete time signal) x=[5 6 47 7 6 -29 4], where x(1)=s(0), x(2)=s(T), x(3)=s(2T),
and so on. The third and sixth elements of x appear to be unusual, and they are called outliers.
Let y denote the number sequence produced by median filtering x. For element n of y, use x(n-
1), x(n) and x(n+1). This is called a K=3 element window, w = [x(n-1) x(n) x(n+1)]. Sort w into
ascending values to obtain a K=3 element vector u. Now set y(n) to the middle element of u.
Then, increment n by 1, get the next w, and repeat to get the next y. For example,
n=2, w=[5 6 47] → sort → u=[5 6 47], y(2)=6
n=3, w=[6 47 7] → sort → u=[6 7 47], y(3)=7
n=4, w=[47 7 6] → sort → u=[6 7 47], y(4)=7
n=5, w=[7 6 -29] → sort → u=[-29 6 7], y(5)=6
n=6, w=[6 -29 4] → sort → u=[-29 4 6], y(6)=4
n=7, w=[-29 4 4] → sort → u=[-29 4 4], y(7)=4. There is no x(8); x(7) was repeated.
the signal x. Median filtering works best to replace outliers with a neighboring element when
other consecutive elements in x have similar values.
We could have used K=4, in which case w is set to w=[x(n-2) x(n-1) x(n) x(n+1)], and u will not
have a middle element. In this case, set y(n)=(u(n-1)+u(n))/2. Or, we can apply the K=3 median
filter again to y, and the result will be even less noisy and most likely too different from x
without the noise. The built-in MATLAB function median finds the median of a vector w.
highest frequency tone. The frequency (pitch) f of the tone produced by key n is given by
the program uses (4.8) to get the frequency of the note. A sinusoid at this frequency is then
sampled at the rate of 44,100 samples/sec for 5 secs. To make the sound produced more
realistic, the sampled sinusoid is multiplied by the pulse shown in Fig. 3.3, which becomes the
envelope of the sound. Therefore, the sound increases quickly and decays slowly. To hear the
sound, a two column (two channel stereo) matrix must be set up, and then, the built-in
MATLAB functions sound or soundsc can be used to send the matrix to the audio play device of
the computer. You can also use the built-in MATLAB function wavwrite, which creates a
Windows WAV file that can be played by any multimedia application. The syntax for wavwrite
is wavwrite(x,fs,’name.wav’), where x is a two column matrix, fs is the sampling rate (also called
sampling frequency), and name can be replaced by a name you can choose.
% Keys are identified by an integer n, n = 1 to 88 from left to right
% Only the tone at the fundamental frequency (pitch) will be sounded
clc; clear all
fs = 44100; % this is the sampling rate used to produce audio CDs.
disp(‘Play a note on a virtual piano.’)
n = input(‘enter an integer from 1 to 88 of a piano key: ‘); % using n=40
n = round(n); % make sure entry is an integer
if n < 1, n = 1; end % keep n in range
if n > 88, n = 88; end
f = 440*(2^(1/12))^(n-49); % the frequency in Hz of the key
w0 = 2*pi*f; % convert Hz to rad/sec
T0 = 5; % the tone will be sounded for 5 secs
delta_t = 1/fs; % the time increment between sound sample
N = round(T0/delta_t); % integer number of samples
time = delta_t*[0:N-1]; % N time points where sound is sampled
p_key = sin(w0*time); % get samples of sound, limit amplitude to +/-1
% get an envelope
load(‘pulse.mat’,’x’) % retrieve the pulse that was built in Example 3.15
p_key = x .* p_key; % multiply the sound by an envelope, this is optional
stereo_sound = [p_key’ p_key’]; % form left and right channels
sound(stereo_sound,fs) % send sound to audio play device
will likely contain some spike like noise. The second part of Prog. 4.7 simulates the addition of
to reduce the spike noise contribution to the sound.
noise = 0.1*randn(1,N); % using the normal random number generator
p_key = p_key + noise; % add noise to signal
% generating randomly positioned spikes
increment = 100; % add noise spike somewhere within every 100 sample of p_key
spikes = zeros(1,N); % zero out and preallocate space
for k = 1:increment:N;
location = floor(k + rand*increment); % use rand to locate spike
if rand > 0.5; % use rand to determine sign of the noise spike
sign = +1.0;
else
sign = -1.0;
end
spikes(location) = sign*2.0*rand; % use rand to determine spike value
end
p_key = p_key + spikes; % add spike noise to signal
plot_points = 1500; range = 5e4:5e4+plot_points-1; % plotting few points
t_interval = time(range); y = p_key(range);
plot(t_interval,y);
axis([t_interval(1) t_interval(plot_points) -3 3])
title(‘Segment of note middle C plus random noise plus random spikes’);
xlabel(‘seconds’);
ylabel(‘amplitude’);
grid on
max_value = max(abs(p_key));
% limit sound to +/- 1, required by audio device
normalized_p_key = p_key/max_value;
stereo_sound =[normalized_p_key’, normalized_p_key’];
sound(stereo_sound,fs)
% median filtering sound
filtered_p_key = zeros(1,N); % preallocate space
K = 3; % use 3 element median window
for k = 1:N
if k == 1
w = [p_key(1) p_key(1) p_key(2)];
filtered_p_key(1) = median(w);
elseif k == N
w = [p_key(N-1) p_key(N) p_key(N)];
filtered_p_key(N) = median(w);
else
w = [p_key(k-1) p_key(k) p_key(k+1)];
filtered_p_key(k) = median(w);
end
end
y = filtered_p_key(range); % plot a segment of filtered_p_key sound
figure % causes function plot to open a new figure window
plot(t_interval,y)
axis([t_interval(1) t_interval(plot_points) -3 3])
title(‘Segment of median filtered note middle C’);
xlabel(‘seconds’);
ylabel(‘amplitude’);
grid on
max_value = max(abs(filtered_p_key));
normalized_filtered_p_key = filtered_p_key/max_value;
stereo_sound =[normalized_filtered_p_key’, normalized_filtered_p_key’];
sound(stereo_sound,fs)
signal is awful. Fig. 4.7 shows the result of median filtering to remove spike noise.
-3
0
1
2
Segment of note middle C plus random noise plus random spikes
pl
de
-3
-2
-1
0
1
2
Segment of median filtered note middle C
am
pl
itu
de
processing can be performed to reduce the noise further.
__________________________________________________________________
expression is no longer true. Its syntax is
include additional while loops and for loops. While loops can be nested. The block of
statements can include if-elseif-else-end structures and any other valid MATLAB statements. It
is useful to indent the block of statements to easily see the matched while and end key words.
The expression must evaluate to a logical or real scalar. If it evaluates to a real number, then it
is treated as a logical value.
2) if the expression is true, then execute the block of statements
3) and go back to the while statement and continue with step (1)
4) if the expression is not true, then skip the block of statements, and continue execution after
the end statement.
while n <= N
n = n +1;
end
Here, the statements within the while loop will be executed N times. However, this is not the
intention of including a while loop structure in addition to a for loop structure in the
programming language. Usually, a while loop structure is used when we do not know ahead of
time how many loop iterations are necessary to accomplish some task. For example,
err = 0;
while ~err
block of statements
unchanged and the loop continues to iterate or set err to err=1, which causes the while loop to
terminate.
while ~negligible
sum_T = sum_T + next_term; % accumulate terms
% prepare the next term in the Taylor series
factorial_term = factorial_term*(2*k)*(2*k+1); % update (2k+1)!
% uses only two multiplies to find each (2k+1)!
minus_one_term = – minus_one_term; % update (-1)^k
x_power_term = x_power_term*x_square; % update x^(2k+1)
% uses only one multiply to find each x^(2k+1)
next_term = minus_one_term*x_power_term/factorial_term;
% compare next term in Taylor series to sum and check precision
if abs(next_term) < abs(precision*sum_T);
sum_T = sum_T + next_term; % include the next term
negligible = 1; % no need for additional loop iterations
end
k = k + 1; % count iterations required to achieve precision
end
The if-end block of statements checks a condition that can change the result of evaluating the
while expression.
______________________________________________________________
The break and continue commands can also be used in while loops. They work the
same way as in for loops. If a while loop expression is always true, then, for example, we have
while 1
block of statements
end
The block of statements must include a test of some condition that can cause a break,
otherwise the loop will execute endlessly. You can terminate a MATLAB program by activating
the Command Window and depressing Ctrl and C on the keyboard.
4.6 Method of Steepest Descent
( )f x has a minimum. Denote the value of x that minimizes ( )f x by optx . Therefore,
( ) ( )optf x f x≤ , for optx x≠ . A method to find optx is based on the drawing shown in Fig. 4.8.
optimizer of ( )f x ), and the other two minima at ax x= and bx x= are called local minima.
at the guess nx is positive, then we know that we must decrease this guess. Similarly, if nx is to
the left (but not too far) of optx , then the slope of ( )f x at nx will be negative, and we know
that we must increase this guess. Therefore, given a guess nx of optx , adjust the value of nx to
obtain the next guess 1nx + with
derivative of ( )f x , assuming it exists, evaluated at nx x= .
where 0x = initial guess of optx . If 0x is within the vicinity of optx and µ is small enough, then
lim n optn x x→∞ = , since as n optx x→ , ‘( ) 0nf x → . Selecting a suitable value of µ is a matter of
step size is too big, then the next value 1nx + may be further away from optx than nx , or even in
converge to optx . If the step size is too small, then (4.9) will require many iterations (much
3) if ε− ≥new oldx x , then old newx x= and repeat step (2), else opt newx x≅ and stop.
of steepest descent, or another kind of optimization method can be employed, if we covert the
given problem to an optimization problem. If it exists, we can find the solution, call it optt , by
2′( ) 2( )( 2 )t tf t e t e t− −= − − − , Prog. 4.8 finds optt . The results follow the program.
% algebraic equation: exp(-t) = t^2
clear all; clc; format long e
precision = 1e-8; % desired accuracy
mu = 1e-4; % step size
t_old = 1.0; % initial guess of the optimizer
f_of_t_initial_guess = (exp(-t_old)-t_old^2)^2
f_prime = 2*(exp(-t_old)-t_old^2)*(-exp(-t_old)-2*t_old);
t_new = t_old – mu*f_prime;
n_iterations = 1; % initialize an iterations counter
max_iter = 1e5; err = 0; % limit the number of iterations
while abs(t_new – t_old) >= precision % check for convergence
n_iterations = n_iterations + 1; % count iterations of algorithm
if n_iterations > max_iter, err = 1; break; end
t_old = t_new; % update recursion
% get new solution
f_prime = 2*(exp(-t_old)-t_old^2)*(-exp(-t_old)-2*t_old);
t_new = t_old – mu*f_prime;
end
if err == 0
t_opt = t_new
f_of_t_opt = (exp(-t_opt)-t_opt^2)^2
fprintf(‘Algorithm required: %6i iterations to converge\n’,n_iterations)
else
fprintf(‘Algorithm exceeded: %6i iterations\n’,max_iter)
t_opt = 7.034812353615486e-001
f_of_t_opt = 6.900847459766987e-010
Algorithm required: 13347 iterations to converge
_____________________________________________________________
be improved by increasing the step size. It requires a derivative of the function to be
optimized, which for some problems may be practically impossible to obtain. However, it and
its variations are widely applied, because it is a simple algorithm to implement. The required
derivative could be approximated with, for example
f x h f x h
h
function evaluations, which can increase computing time. Also, once nx is close to the
optimizer, (4.11) will give small and erratic ‘( )nf x values, making further convergence unlikely.
However, this depends on the behavior of ( )f x and the desired precision.
free optimization method [1], which makes it attractive to use. A syntax is given by
of a function m-file, or it could be a handle of an anonymous function. Use doc fminsearch for
more syntax options.
given to define ( )f t , as follows.
>> f_of_t = @(t)(exp(-t)-t^2)^2; % define an anonymous function
>> f_of_t(1)
ans = 3.995764008937280e-001
>> [t_opt f_opt] = fminsearch(f_of_t,1)
t_opt = 7.035156249999998e-001
f_opt = 8.403998856470625e-009
function evaluations or never converge to useful results. This depends on the behavior of the
function to be optimized and the initial guess. The function ( )f t of Example 4.14 has only one
minimum, the global minimum, and both methods of optimization used here worked very well,
regardless of the initial guess. However, there are no known methods of optimization that can
with certainty find the global minimum of a function. This depends on the behavior of the
function, and it is important to provide an initial guess in the vicinity of the optimizer.
t a t
−∞ −∞
t
d q t
d t
t t a t
−∞ −∞
from t =− ∞ until t a= , and the first term is the additional amount of charge that has passed
through the wire cross-section from t a= until sometime t , where t a> .
a set of discrete time points.
Consider the integral
T
interval is 0 /T T N= . This is depicted in Fig. 4.9. The value of ( )x t at the beginning of each
sub-interval is ( ), 0, 1, , 1x nT n N= − . Euler’s method approximates the area under the
curve in each sub-interval by a rectangular area given by the width T times the value of ( )x t at
the beginning of the sub-interval. The integral of (4.15) is approximated by summing all of the
rectangular sub-interval areas to get
A x kT T T x kT
approximation gets better as N is increased, which makes T smaller.
( 1)N T−
( )x T
(( 1) )x N T−
( ) ( )
S n x k T
=
and then ( )A T S N . For some n , assume that we have ( ), 0, , 1x kT k n= − , which means
that we can find ( )S n . When ( )x nT becomes available, we can update ( )S n to get ( 1)S n +
with
S(2) = S(1) + ( )
recursion (4.17) is operated for 0, 1, , ( 1)n N= − , with (0) 0S = .
, 0
x t t
t
≠
With
>> x = sinc(t); % evaluate x(t) = sin(π t)/π t at the discrete time points given by the elements of t
>> plot(t,x); grid on
>> xlabel(‘t – secs’); ylabel(‘sinc(t)’);
with an amplitude that decays to zero. The integral of ( ) sin c( )x t t= from −∞ to +∞ is one.
-0.5
0.5
1
nc
)
0.4898881. Therefore, the integral from 10t = to t →+ ∞ is (0.5 – 0.4898881). Prog. 4.9 uses
Euler’s method to find this integral, and the results for N = 50, 550, 1050 follow.
clear all; clc
T0 = 10; % time range is from 0 to T0
for m=1:3
N = 50+500*(m-1); % number of sub-intervals
T = T0/N; % sub-interval width
S = 0; % initialize the sum
for n = 0:N-1 % loop to include each sub-interval
t = n*T; % time point
S = S + sinc(t); % recursion
end
Area(m) = T*S; % area
N_subs(m) = N;
end
format short
table = [N_subs;Area];
disp(‘ N Area’)
fprintf(‘ %i %f \n’,table)
50 0.590224
550 0.498982
1050 0.494651
Nevertheless, depending on the behavior of ( )x t , Euler’s method is often used because of its
simplicity.
of using the value of ( )x t at the beginning of each sub-interval, the trapezoidal rule uses the
average of the values of ( )x t at the beginning and end of each sub-interval, and (4.16) and
(4.17) become
0 0
( ) (( 1) )
k k
A T x kT x k T
= =
= + +
loop in Prog. 4.9 is replaced by
for n = 0:N-1 % loop to include each sub-interval
t = (n+1)*T; % next time point
% evaluating the sinc function just once for each loop iteration
xn1T = sinc(t);
S = S + xnT + xn1T; % trapezoidal rule recursion
xnT = xn1T; % update x(nT)
end
50 0.490224
550 0.489891
1050 0.489889
______________________________________________________________
these functions is quad. It is a function function, and the syntax to use this function is
my_function.m. The fourth argument specifies an error tolerance. If tol is not included in the
input argument list, then MATLAB uses a default value of 1e-6.
function in Example 4.15. The statement that was used is: area=quad(@sinc,0,10.0,1e-8).
defined function a vector input. Therefore, the user-defined function must be properly
vectorized to return a vector. See Table 3.3 for more details about vectorizing expressions
within a function that involve function input vectors.
sinusoidal voltage across it. Recall the diode i-v characteristic, for which a function is given
below.
% current through a diode when the voltage across it is v=Asin(wt)
global w A % frequency and amplitude
v = A*sin(w*t);
I_S = 1e-12; % saturation current
V_T = 25.85e-3; % thermal voltage at 300 K
current = I_S*(exp(v/V_T) – 1);
sinusoidal voltage. The frequency and amplitude of the voltage are passed to this function as
global variables. This function was saved as diode.m. To see how charge moves through the
diode over one cycle of the voltage, let us segment one cycle into N sub-intervals, and
integrate the current over each sub-interval. Prog. 4.10 gives a script to plot charge movement
over one cycle.
% a sinusoidal voltage across it
clear all; clc
global w A
w = 2*pi*frequency;
A = 0.8; % amplitude of sinusoid
T0 = 1/frequency; % time of one cycle
N = 100; T = T0/N; % will use N+1 plot points
tol = 1e-8; % specify tolerance
q(1) = 0; % initial charge is zero
v(1) = 0; % initial voltage value
t = 0.0:T:T0; % specify N+1 time points
for k = 1:N
v(k+1) = A*sin(w*t(k+1));
q(k+1) = q(k) + quad(@diode,t(k),t(k+1),tol); % accumulate charge
end
plotyy(t,q,t,v); grid on % use right y-axis for v
xlabel(‘time secs’); ylabel(‘Coulombs’)
title(‘Charge movement through a diode’)
0
ulo
bs
–
0.2
0.4
0.6
0.8
it. Notice that as the voltage increases, there is no appreciable movement of charge through
the diode until the voltage reaches a threshold voltage of about 0.7 volts. Then, once the
voltage exceeds this threshold, the amount of charge that has moved through the diode
increases quickly until the voltage again goes below the threshold, after which the net amount
of charge through the diode no longer increases.
______________________________________________________________
adaptively adjust the way they work to achieve prescribed error tolerances. For example, quadl
smoothly behaving integrands. The built-in function quad2d numerically obtains double
integrals over a planar surface. See the MATLAB help facility for details and additional built-in
functions for numerical integration.
set of blocks of statements should be executed. The syntax of this structure is
case expression_1
#N of statements is optional. The switch expression must evaluate to a scalar or a character
string. When MATLAB encounters the switch key word, it does the following
2) if the result of evaluating the switch expression matches the case expression_1 value, then
execute block #1 statements, and continue execution after the end statement
3) if the result of evaluating the switch expression does not match the case expression_1 value,
then skip the block #1 statements to the next case
4) if the result of evaluating the switch expression matches the case expression_2 value, then
execute block #2 statements, and continue execution after the end statement
5) if the result of evaluating the switch expression does not match the case expression_2 value,
then skip the block #2 statements to the next case, and so on
6) if the result of evaluating the switch expression does not match any case expression value,
then execute the block #N statements after the otherwise key word (if there is one) and
continue execution after the end statement.
that each switch key word is matched by an end key word. Within each block of statements
there can be any of the script flow control structures, including nested switch-case-otherwise
structures.
depending on the color and behavior of the traffic light. Prog. 4.11 is a demonstration of the
application of the switch-case-otherwise structure to this situation.
% Local conditions may require other actions
clear all; clc
disp(‘Upon approaching an intersection with working traffic lights,’)
disp(‘different actions are required.’)
behavior = input(‘Enter a traffic light behavior in single quotes: ‘);
switch behavior
case ‘red’
disp(‘Come to a complete stop; and wait for green light to proceed.’)
case ‘yellow’
disp(‘Slow down; prepare to stop.’)
case ‘green’
disp(‘Proceed with caution.’)
case ‘flashing yellow’
disp(‘Slow down; proceed with caution; be ready to stop.’)
case ‘flashing red’
disp(‘Come to a complete stop; proceed with caution.’)
otherwise
disp(‘Traffic lights may be out of service.’)
disp(‘Come to a complete stop.’)
disp(‘Local conditions may require other actions.’)
end
script are essential for implementing problem solution methods. All programming languages
have means to do this. You should now know how:
• to use logical operators to make compounded statements using logical parts
• to make a script more robust
• the for loop structure works to repeatedly execute a block of MATLAB statements
• to find the time required to execute a block of statements
• a histogram is obtained
• to determine the probability of an event by counting
• a pseudorandom number generator can be used to simulate random events
• a median filter works
• to use the sound device
• the while loop structure works to repeatedly execute a block of statements depending
• Euler’s method and the trapezoidal rule are used for numerical integration
• the switch-case-otherwise structure is used to execute alternative script parts
MATLAB’s help facility to learn more about these built-in functions, where you will also find out
about many other related built-in functions.
isequal(A,B) returns logical 1, meaning true, if arrays A and B are the same size
otherwise
contain a non-zero element at that same array location, and
otherwise, that element is set to 0
contains a non-zero element at that same array location, and
otherwise, that element is set to 0
contains a non-zero element at that same array location, but not
both, and 0 if both elements are zero or non-zero
value
element at that same array location, and 0 otherwise
|| short logical or
I = find(A) returns the linear indices corresponding to the non-zero entries of
L = sign(A) signum function, returns 1 if the element is greater than zero, 0 if
elements of complex A, sign(A) = A ./ abs(A)
tic
fraction; just pause causes a script or function to stop and wait for
the user to strike any key before continuing
bars; bar(Y) uses the default value of X=1:M. For vector inputs,
bar(X,Y) or bar(Y) draws length(Y) bars. There is no space
between bars if width=1.
break terminate execution of for or while loop
return causes a return to the invoking script or function or to the
keyboard when placed in an m-file, stops execution of the file and gives
K appearing before the prompt. Variables may be examined or
changed; all MATLAB commands are valid. The keyboard mode is
terminated by executing the command return.
uniformly distributed over the open interval (0,1)
uniformly distributed over the range 1:imax
the standard normal distribution
Y = var(X) Y is a row vector containing the variance of each column of X
I = hist(X,N) bins the elements of X into N bins
fprintf write formatted data to a display or a file; more about this later
y = median(X) returns the median value of the elements in the vector X
sound(X,fs) sends the signal in vector X, with sampling frequency fs to the
values outside the range -1 to +1.
wavwrite(X,fs,’file.wav’) writes data X to a Windows WAV file specified by the file name
as a matrix with two columns
elements in the matrix X
function specified by the function handle func
versus x2 with y-axis labeling on the right, using the plotting
function specified by the function handle plot_func
method
f(x,y) specified by the function handle func, where ≤ ≤a x b and
four chapters of this book. Use demo to watch an excellent audio/video tutorial about if-elseif-
else, for and while loop control structures. In these chapters some special features, for
example, plotting, were introduced to make examples more complete with visual displays. In
the remaining chapters, special, unique and particularly useful features of this programming
language will be discussed.
role in Boolean algebra, logic circuits and binary arithmetic.
7, pp. 308–313.
J. for Numerical Methods in Engineering, vol. 60, pp. 1301-1315.
state the purpose and activity of the script.
d = a^2 >= c, (b) e = (a > b -4)+c <= b^2, (c) f = a+b == c-4.
2) For t = -1:0.1:2, w = 2*pi and x = cos(w*t), use a logical operation to find a vector y, which has
the same dimension as x, where all of the positive elements of x have been replaced by zero.
3) A = [-1 9 2; 0 -3 4], B = [], C = ones(3,2)'. Give the result of: (a) isempty(A), (b) isempty(B), (c)
logical(A), (d) any(A), (e) islogical(C), (f) C < false(2,3) | A, (g) isinteger(sqrt(A)), (h) isfloat(C).
Section 4.2
4) For a = 2, b = 0 and c = -5, give the results of: (a) d = a | (b & c), (b) e = and(a,b)|c,
(c) f = a & c | ( b & ~a), (d) g = a & c | b & ~a, (e) h = (b && c) | (b || a), (f) p = a & xor(~b,c).
5) (a) Describe the difference between (a & b) and (a && b).
(b) Describe the difference between (a | b) and (a || b).
(b) e = sqrt(a^2 + b^2) < abs(c) | (a+b)/2 >= c ,
(c) f = ~isempty(find([a b 2*c sqrt(a^2 + c^2)] > 4)).
y. (b) let B = reshape(A,1,9) and give C = B(y).
(b) If the inputs are the same as in part (a), except that x = 2.7, then what is the script output?
b, x and N, and the output arguments are: k and err. Include comments that explain the
purpose of the function and the meaning of the err codes. Provide a copy of the function.
(b) Write a script that demonstrates using the function bin_number. Provide a script listing
and demonstration of script operation.
structure. Provide a script listing and a demonstration of script operation. Does the current
have to be a positive number?
Display message
NoYes
delivered to resistor
power
input is an integer and that it is within the required range. Otherwise, display an error message
for each input error type. Then, the script must determine if N is an odd or even integer and
display a message about the result. Provide a script listing and a demonstration of script
operation.
what will MATLAB do when the for loop is completed with an end statement?
A_dot_B of two vectors A and B with dimension 1xN. The function input arguments must be A
and B, and the function returns A_dot_B.
(b) Write a script to demonstrate the operation of your function where t = 0:T:T0-T, T = T0/N
secs, T0 = 02 /π ω secs, A = cos( 0ω t) and B = cos(K 0ω t) for K = 1, 2 and 5. You choose N and
(d) When the inner product of two vectors is zero, it is said that the vectors are orthogonal.
Here, we find that this concept can be extended to functions, e.g., cos(M 0ω t) and cos(K 0ω t),
M ≠ K, are orthogonal functions over a period T0. What other sinusoidal functions are
orthogonal over a period?
series computation when the next term in the series contributes less than 0.0001 % of the sum
of terms. For example, see Prog. 4.4. Provide a figure of a plot of your result. Include in the
figure a stem plot that uses the MATLAB built-in function exp. The figure should include a grid,
axes labels and a title.
(a) In terms of V, Rs and R, obtain an expression for the power P delivered to R.
(b) Write a script that receives as input the values of V and Rs, and makes a vector R with
element values that range from 0.1*Rs to 10*Rs in increments of 0.1*Rs. Use a for loop to
create a vector of power delivered to each element value of R. Then, plot the power P versus R.
Do this for V = 10 volts, Rs = 10 Ohms and Rs = 100 Ohms. For each case, from inspection of the
plots, give a value of R such that V delivers the maximum power to R. Provide a script listing
and plots, including axes labels and titles.
(c) Vectorize your expression for P, and repeat part(b) without using a for loop.
(d) In general, for a given value of Rs, what value of R enables the resistor R to receive the most
power from the voltage source V?
numbers with the function randn.
(a) Obtain the mean, variance and histogram of x. Use 51 bins for the histogram. Use the
MATLAB built-in functions mean, var and hist to obtain these results.
(b) Let y = ax +b, where a = 3 and b = 2. Use a for loop to find the mean of y. Then, use another
for loop to implement (4.6). How does your result compare with using the MATLAB function var
to find the variance of y?
(c) Obtain a histogram of y, using 51 bins. How are the parameters a and b related to the
variance and mean of y?
6. This can be simulated with die_1 = ceil(6*rand) and die_2 = ceil(6*rand), where each
invocation of rand produces a random number from 0 to 1. Write a script that uses a for loop to
do N_exps = 1e7 experiments, and within the loop uses an if-elseif-end structure to count the
number of times die_1 + die_2 = 7 or 11. Begin your script with the statement: rng(‘default’) to
bring the random number generator to its default initial state. By counting, find the probability
that an experiment outcome will be 7 or 11.
4
1 2( ) sin( ) 0.5 cos( )s t t tω ω= + that is corrupted by additive noise ( )n t . Let 1 2 (262)ω π=
data ( )x k T for this simulation is given by ( ) ( ) ( )x kT s kT a n kT= + , where ( )n kT is found
by using a normal random number generator and a controls the noise amplitude. To help you
get started, use
% Program to implement a digital low-pass filter
clear all; clc;
fs =8000; T = 1/fs;
w1 = 2*pi*262; w2 = 2*pi*440;
a = 0.5; % used to see how noise affects the result
k = 0; t = k*T; time(k+1) = t; % initialize sample time
x(k+1) = sin(w1*t) + 0.5*cos(w2*t) + a*randn; % first sample
% since there is no x(-T) and x(-2T) use
y(k+1) = x(k+1)/4;
k = k+1; t = k*T; time(k+1) = t; % save time for plotting
x(k+1) = sin(w1*t) + 0.5*cos(w2*t) + a*randn; % second sample
y(k+1) = (x(k+1)+2*x(k))/4;
N = 100; % process N-2 more samples of the input
for k = 2:N-1
t = k*T; time(k+1) = t; % save time for plotting
x(k+1) = sin(w1*t) + 0.5*cos(w2*t) + a*randn; % kth sample
y(k+1) =
Complete the for loop with the given algorithm, and then plot x and y versus time. What does
the algorithm appear to do?
until a program user enters an R value less than or equal to zero. For example,
while 1
R = input(‘Enter a resistance value: ‘)
if R <= 0
disp('Terminating program')
break
end
Complete the script. For a nice looking output, use, for example, fprintf('Power delivered =
%6.2f watts\n', P). Provide a script listing, and demonstrate its operation.
21) Do Prob. 4.15 using a while loop inside a for loop. With the for loop index ranging from n
= 1 to n = N, where N = length(t), then, for example, the while loop could be
x = -2.0*t(n);
exp_sum = 1.0; % first term in power series
next_term = x; k = 1; % second term in power series
while abs(next_term) > 1e-6*abs(exp_sum)
exp_sum = exp_sum + next_term;
k = k+1;
next_term = next_term*x/k;
end
exp_func(n) = exp_sum + next_term; % use last term
described in Prob. 4.15.
find the value of R that maximizes the power delivered to it by the voltage source V. Use input
statements to obtain the values of V and Rs, and a while loop to implement an update method
to find Rnew given Rold similar to the update method of (4.9). See Example 4.14. Find a suitable
value of the step size by trial and error. Use Rold = Rs/2 as the initial guess of Ropt. Provide a
script listing and a demonstration of its operation.
designed to find y. This can be done by converting the given problem to an optimization
problem. Consider minimizing the function f(y) = (y2 – x)2, where f ‘(y) = 2(y2-x)(2y). Write a
script that uses a while loop within a while loop to apply the method of steepest descent. The
outer while loop must contain an input statement to receive a value of x, and if the input is not
greater than zero, then break out of this loop, which terminates the script. If the input x is
greater than zero, then the inner while loop must apply the method of steepest descent to find
y. Implement an update method to find ynew given yold like the update method of (4.9). Use yold
= x as an initial guess of yopt. See Example 4.14. Find a suitable step size by trial and error.
Provide a script listing and a demonstration of its operation for three positive values of x.
function fminsearch. Create a function m-file to define the function to be optimized by
fminsearch. Provide listings of your script, the function m-file, and a demonstration of script
operation.
(a) Manually, find an expression for the area A under x(t) from t = 0 to t = 5.
Write a script to:
(b) plot x(t) for t = 0 to t = 5. Use N = 501 values of t in the given time range
(c) calculate A using the expression found in part(a)
(d) apply Euler’s method to find an approximation of A for N = 11, 111 and 211
(e) output a table of the difference between A and its approximation for each N.
Use format long e for all script output. Provide a listing of your script and the results.
under x(t) from t = 0 to t = 5. Use an anonymous function to define x(t). Compare the result to
your script and results.
column vector Y and a character string named solution_method, and the outputs are the
condition number K of A and the solution X of AX=Y, assuming that the determinant of A is not
zero. First, check the determinant of A, and if it is zero, set solution_method = ‘none’. The
possible solution methods you should use are: inv (the function), left matrix divide, eigenvalues
and eigenvectors or singular value decomposition, as designated with solution_method = ‘inv’,
‘left_mat_div’, ‘eigen’ or ‘singular_decomp’. Use the switch-case-otherwise-end structure to
select the method to be used. Use the otherwise option to return an empty vector X and K = inf
in the event that the character string solution_method is not one of the allowed solution
method options. Write a script that demonstrates the application of your function. Provide
script and function listings and script output.
number of times that an experiment outcome is 7 or 11. Use: die_1 + die_2, for the switch
expression and case 7 and case 11.
