funplot.m

很多matlab的源代码

function [t0,y0] = funplot(fname,lims,npts,angl,subdiv)
% FUNPLOT Plot of string function of t.
%
%	FUNPLOT(FNAME,L) plots FNAME over limits in L = [tmin tmax].
%	FNAME may be a function m-file such as 'sinc' OR 
%	a string function of t such as 'sin(t)' or '2*t.*exp(-3*t)'.
%
%       FUNPLOT(FNAME,L,N) uses a minimum of N samples [Default: 101]
%
%       FUNPLOT(FNAME,L,N,D) plots more samples if contiguous segments are
%       at an angle larger than D degrees [Default: 5].
%
%       FUNPLOT(FNAME,L,N,D,S) uses no more than S iterations for rapidly
%       changing portions [Default: 7].
%
%       [T,Y] = FUNPLOT(...) returns T & Y data, no plot is drawn.
%
%       USAGE: An example using defaults: funplot('exp(-t).*sin(5*t)',[0 4])
%
%       FUNPLOT (with no input arguments) invokes the following example:
%
%	% Plot sinc(t)=sin(pi*t)/(pi*t) over [-4 4].
%         >>funplot('sinc',[-4 4])  %Note: sinc.m is a function m-file 

% Modified version of FPLOT to accomodate string functions
% Modified by A.Ambardar 01-04-93. With permission from The Mathworks, Inc.


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998


if nargin==0,help funplot,disp('Strike a key to see results of the example')
pause,funplot('sinc',[-4 4]);return,end

if nargin < 2
    error('Plotting limits [lo hi] must be specified.');
elseif nargin < 5
    subdiv = 7;
    if nargin < 4
        angl = 5;
    end
    if nargin < 3
        npts = 101;
    end
end
if isstr(lims) | isstr(npts) | isstr(angl) | isstr(subdiv)
        error('Optional arguments must be numeric.');
end
iter = 0;
angl = angl*pi/180;
angl = abs(angl);
ang = angl + 1;
tmin = min(lims);
tmax = max(lims);
t = tmin + (0:npts-1)'*(tmax-tmin)/(npts-1);

%%%NEW ADDITION & CHANGE Ambardar
ff=0;if exist(fname)>1,y=feval(fname,t);else,ff=1;y=eval(fname);end

ynan = y(~isnan(y));
yscal = max(max(abs(ynan)));
y = y/yscal;
close_enough = 0;

[m,n] = size(y);
isvector = (n~=1);

% see if need to fill in values for 'smoothness'
while any(ang > angl) & (iter <= subdiv) & (~close_enough)
    [m,n] = size(y);
    yi = diff(y);
    ti = diff(t)*ones(1,n);
    nt = length(ti) - 1;
    radi = sqrt(ti.^2+yi.^2);
    num = ti(1:nt,:).*ti(2:nt+1,:) + yi(1:nt,:).*yi(2:nt+1,:);
    den = radi(1:nt,:).*radi(2:nt+1,:);
    ang = abs(acos(num./den));
    if isvector, ang = max(ang.').'; end
    iter = iter + 1;
    ii = find(ang > angl);    % need to interpolate more
    if isempty(ii)
        break;
    else
        % Choose new points 1/3 to each side of points with large angle.
        tt = [(t(ii) + 2*t(ii+1))/3; (2*t(ii+1) + t(ii+2))/3]; 


%%%NEW ADDITION & CHANGE Ambardar
if ff==0,yt=feval(fname,tt)/yscal;
else,t1=t;t=tt;yt=eval(fname)/yscal;tt=t;t=t1;end

if all(all((abs([y(ii,:);y(ii+1,:)]-yt)<.01) & ...
(abs([y(ii+1,:);y(ii+2,:)]-yt)<.01)))
                close_enough = 1;
        end
        [t, ind] = sort([t;tt]);
        y = [y;yt];
        y(:) = y(ind,:);
    end
end
y = y*yscal;
if iter >= subdiv
    disp('WARNING: Iteration limit reached.');
    if nargout == 0
        disp('         Plot may be inaccurate.');
    end
end
if nargout == 0,
if exist('version')==5,eval('animate(t,y)');else,plot(t,y,'r'),grid,end
else
    t0 = t;
    y0 = y;
end