Applied Numerical Methods Using MATLAB. Won Y. Yang
or ‘ fplot()’? Do they work properly?>ezplot('tan(x)',[0 2*pi]) >fplot(@(x)tan(x),[0 2*pi]) % plot an anonymous functionFigure P1.5 Plotting the graph of f(x) = tan x.
6 1.6 Plotting the Graph of a Sinc FunctionThe sinc function is defined as(P1.6.1) whose value at x = 0 is(P1.6.2) We are going to plot the graph of this function over [ ‐4π,4π].(a) Casually, you may try as follows:>x=[-100:100]*pi/25; y=sin(x)./x; >plot(x,y), axis([-15 15 -0.4 1.2]) Even if no warning message shows up, is there anything odd about the graph?(b) How about trying with a different domain vector?>x=[-4*pi:0.1:+4*pi]; y=sin(x)./x; >plot(x,y), axis([-15 15 -0.4 1.2]) Surprisingly, MATLAB gives us the function values without any complaint and presents a nice graph of the sinc function. What is the difference between (a) and (b)?(cf) Actually, we would have no problem if we used the MATLAB built‐in function ‘ sinc()’.
7 1.7 Term‐Wise (Element‐by‐Element) Operation in MATLAB User FunctionsLet the function f1(x) be defined without one or both of the dot ( .) operators in Section 1.1.6. Could we still get the output vector consisting of the function values for the several values in the input vector? You can run the following statements and see the results.>f1=@(x)1./(1+8*x̂2); f1([0 1]) >f1=@(x)1/(1+8*x.̂2); f1([0 1]) Let the function f1(x) be defined with both of the dot ( .) operators as in Section 1.1.6. What would we get by running the following statements?>f1=@(x)1./(1+8*x.̂2); f1([0 1]')
8 1.8 In‐line Function and M‐file Function with the Integral Function ‘ quad()’As will be seen in Section 5.8, one of the MATLAB built‐in functions for computing numerical integrals is ‘ quad()’, the usual usage of which is(P1.8.1) wheref: the name of the integrand function (M‐file name should be categorized by ' ' )a,b: the lower/upper bound of the integration intervaltol: the error tolerance (10−6 by default []), trace: 1(on)/0(off) (0 by default [])p1,p2, ..: additional parameters to be passed directly to function fLet us use this ‘ quad()’ function with an in‐line function and an M‐file function to obtain(P1.8.2a) and(P1.8.2b) where(P1.8.3) Complete and run the following script “nm01p08.m” to compute the two integrals (P1.8.2a and P1.8.2b).%nm01p08.m m=1; sigma=2; x0=1; Gpdf='exp(-(x-?).̂2/2/?????̂2)/sqrt(2*pi)/sigma'; xGpdf=inline(['(x-x0).*' Gpdf],'x','m','?????','x0'); x2Gpdf=inline(['(x-x0).̂2.*' Gpdf],'x','?','sigma','x0'); int_xGpdf=quad(xGpdf,m-10,m+10,[],0,?,sigma,x0) int_x2Gpdf=quad(x2Gpdf,m-10,m+10,[],0,m,?????,x0)
9 1.9 μ‐Law Function Defined in an M‐FileThe so‐called μ‐law function and μ−1‐law function used for uniform quantization is defined as(P1.9.1a) (P1.9.1b) Below are the incomplete μ‐law function ‘ mulaw()’, μ−1‐law function ‘ mulaw_inv()’, and a script “nm01p09.m” that are all supposed to be saved as M‐files. Complete them and run the script to do the following jobs with μ = 10, 50, and 255:– Find the values y of the μ‐law function for x=[‐1:0.005:1] and plot the graph of y vs. x.– Find the values x0 of the μ−1‐law function for y.– Compute the discrepancy between x and x0.function [y,xmax]=<b>mulaw</b><![CDATA[(x,mu,ymax) if nargin<3, ymax=1; end xmax=max(abs(x)); y=ymax*log(1+mu*???(x/xmax))./log(1+mu).*????(x); % Eq.(P1.9a)function x=<b>mulaw_inv</b><![CDATA[(y,mu,xmax) % Inverse of mu-law if nargin<3, xmax=1; ymax=1; else ymax=max(abs(y)); end if nargin<2, mu=255; end x=xmax.*(((1+??).̂(abs(?)/y???)-1)/??).*sign(y); % Eq.(P1.9b)%nm01p09.m: to plot the mulaw curve x=[-1:.005:1]; mu=[10 50 255]; for i=1:3 [y,xmax]=mulaw(x,mu(i),1); x0=mulaw_inv(y,mu(i),xmax); discrepancy=norm(x-x0) plot(x,y,'b-', x,x0,'r-'), hold on end
10 1.10 Analog‐digital converter (ADC)Below are two ADC routines ‘ adc1(a,b,c)’ and ‘ adc2(a,b,c)’, which assign the corresponding digital value c(i) to each one of the analog data belonging to the quantization interval [ b(i), b(i+1)]. Let the boundary vector and the centroid (level) vector be, respectively, b=[-3 -2 -1 0 1 2 3]; c=[-2.5 -1.5 -0.5 0.5 1.5 2.5];Note that the digital value corresponding to an analog data smaller than b(1) (the smallest boundary) should be c(1) (the smallest centroid) and that corresponding to an analog data larger than b(N) (the largest boundary) should be c(N) (the largest centroid) where the boundary and centroid vectors are assumed to be arranged in ascending order and N is the number of the centroids (quantization levels) or quantization intervals.Explain the jobs of lines 8, 9, and 10 of the following ADC function ‘ adc1()’.Explain the jobs of lines 3, 5, and 7 of the following ADC function ‘ adc2()’.Run the script “nm01p10.m” (see below box) to get two graphs as shown in Figure P1.10. Explain about the graphs.Figure P1.10 The characteristic of an analog‐to‐digital converter (ADC).function d=<b>adc1</b><![CDATA[(a,b,c) %Analog-to-Digital Converter %Input a=analog signal, b(1:N+1)=boundary vector % c(1:N)=centroid vector %Output: d=digital samples N=length(c); for n=1:length(a) I=find(a(n)<b(2:N)); if ∼isempty(I), d(n)=c(I(1)); else d(n)=c(N); end endfunction d=<b>adc2</b><![CDATA[(a,b,c) N=length(c); d(find(a<b(2)))=c(1); for i=2:N-1 index=find(b(i)<=a&a<b(i+1)); d(index)=c(i); end d(find(b(N)<=a))=c(N);%nm01p10.m b=[-3 -2 -1 0 1 2 3]; % Boundary vector c=[-2.5 -1.5 -0.5 0.5 1.5 2.5]; % Centroid vector % Plot the input-output relationship xa=[-300:300]/100; % Analog data in the range [-3∼+3] xd1=adc1(xa,b,c); xd2=adc2(xa,b,c); subplot(221) plot(xa,xd1, xa,xd2,'r:') % Output of ADC to a sinusoidal input t=[0:200]/100*pi; xa=3*sin(t); xd1=adc1(xa,b,c); %xd2=adc2(xa,b,c); subplot(222) plot(t,xa, t,xd1,'r')
11 1.11 Decimal‐to‐Binary/Octal ConversionsReferring to the decimal‐to‐binary conversion process illustrated in Figure 1.6, complete the following function ‘ dec2bin_my()’ so that it can perform the conversion process. Then, using the function, convert two decimal numbers 3 and 14 to their corresponding binary numbers and compare the results with those obtained using the MATLAB built‐in function ‘ dec2bin()’. Can you modify the function so that the binary numbers corresponding to each of x can be placed columnwise, that is, from top to bottom, rather than from left to right when a vector consisting of multiple decimal numbers is given as the first input argument x?Complete the following function ‘ bin2dec_my()’ so that it can perform the binary‐to‐digital conversion process:function y=<b>dec2bin_my</b><![CDATA[(x) % converts given decimal numbers into binary numbers of N bits y=[]; xmax=max(x); N=1; while xmax>1 xmax=floor(xmax/2); N=N+1; end for n=1:length(x) xn=x(n); for i=?:-1:1; xn_1= floor(xn/?); yn(i)= xn-?*xn_1; xn=xn_1; end y= [y yn]; endfunction y=<b>bin2dec_my</b><![CDATA[(x) [M,N]=size(x); % Number of binary numbers and Number of bits for m=1:M y(m,:)=x(m,:)*?.̂[?-1:-1:0]'; endThen, using the function, convert two binary numbers [0 0 1 1] and [1 1 1 0] to their corresponding digital numbers and compare the results with those obtained using the MATLAB built‐in function ‘ bin2dec()’, by typing the following statements at MATLAB prompt:bin2dec_my([0 0 1 1; 1 1 1 0]) bin2dec(['0011'; '1110'])(c) Modify the decimal‐to‐binary converting function ‘ dec2bin_my()’ into another function ‘ dec2oct_my()’ so that it can perform the decimal‐to‐octa conversion process. Then, using the function, convert two decimal numbers 14 and 3 to their corresponding octal numbers and compare the results with those obtained using the MATLAB built‐in function ‘ deci2oct()’.
12 1.12 Truth Table for a Logical (Boolean) ExpressionNoting that the logical AND, OR, and NOT operators are &, |, and ∼ in MATLAB, complete and run the following script “make_truth_table.m” to construct the truth table for the following logical (Boolean) expression:(P1.12.1) that lists the (output) value of the expression for all possible values of the input variables.%make_truth_table.m % makes the truth table for a logical (Boolean) expression N=3; % Number of Boolean variables Inputs=dec2bin_my([0:2̂N-1]); %All possible values of input variables A=Inputs(:,1); B=Inputs(:,2); C=Inputs(:,3); truth_table=[Inputs (?A&?&C)?(B?∼C)]
13 1.13 Playing with PolynomialsPolynomial evaluation: polyval() Write a MATLAB statement to compute(P1.13.1) (b) Polynomial addition/subtraction by using compatible vector addition/subtractionWrite a MATLAB statement to add the following two polynomial coefficient vectors:(P1.13.2) (c) Polynomial multiplication: conv()Write