雷达地杂波仿真_空时自适应信号处理原理_对雷达杂波进行仿真_matlab源码

function [Q,cnt] = quad_m(funfcn,a,b,tol,trace,varargin) %QUAD Numerically evaluate integral, low order method. % Q = QUAD('F',A,B) approximates the integral of F(X) from A to B to % within a relative error of 1e-3 using an adaptive recursive % Simpson's rule. 'F' is a string containing the name of the % function. Function F must return a vector of output values if given % a vector of input values. Q = Inf is returned if an excessive % recursion level is reached, indicating a possibly singular integral. % % Q = QUAD('F',A,B,TOL) integrates to a relative error of TOL. Use % a two element tolerance, TOL = [rel_tol abs_tol], to specify a % combination of relative and absolute error. % % Q = QUAD('F',A,B,TOL,TRACE) integrates to a relative error of TOL and % for non-zero TRACE traces the function evaluations with a point plot % of the integrand. % % Q = QUAD('F',A,B,TOL,TRACE,P1,P2,...) allows parameters P1, P2, ... % to be passed directly to function F: G = F(X,P1,P2,...). % To use default values for TOL or TRACE, you may pass in the empty % matrix ([]). % % See also QUAD8, DBLQUAD. % C.B. Moler, 3-22-87. % Copyright (c) 1984-98 by The MathWorks, Inc. % $Revision: 5.16 $ $Date: 1997/12/02 15:35:23 $ % [Q,cnt] = quad(F,a,b,tol) also returns a function evaluation count. if nargin < 4, tol = [1.e-3 0]; trace = 0; end if nargin < 5, trace = 0; end c = (a + b)/2; if isempty(tol), tol = [1.e-3 0]; end if length(tol)==1, tol = [tol 0]; end if isempty(trace), trace = 0; end % Top level initialization x = [a b c a:(b-a)/10:b]; y = feval(funfcn,x,varargin{:}); fa = y(1); fb = y(2); fc = y(3); if trace lims = [min(x) max(x) min(y) max(y)]; if ~isreal(y) lims(3) = min(min(real(y)),min(imag(y))); lims(4) = max(max(real(y)),max(imag(y))); end ind = find(~isfinite(lims)); if ~isempty(ind) [mind,nind] = size(ind); lims(ind) = 1.e30*(-ones(mind,nind) .^ rem(ind,2)); end axis(lims); plot([c b],[fc fb],'.') hold on if ~isreal(y) plot([c b],imag([fc fb]),'+') end drawnow end lev = 1; % Adaptive, recursive Simpson's quadrature if ~isreal(y) Q0 = 1e30; else Q0 = inf; end recur_lev_excess = 0; [Q,cnt,recur_lev_excess] = ... quadstp(funfcn,a,b,tol,lev,fa,fc,fb,Q0,trace,recur_lev_excess,varargin{:}); cnt = cnt + 3; if trace hold off axis('auto'); end if (recur_lev_excess > 1) warning(sprintf('Recursion level limit reached %d times.', ... recur_lev_excess )) end %------------------------------------------------------------ function [Q,cnt,recur_lev_excess] = quadstp(FunFcn,a,b,tol,lev,fa,fc,fb,Q0,trace,recur_lev_excess,varargin) %QUADSTP Recursive function used by QUAD. % [Q,cnt] = quadstp(F,a,b,tol,lev,fa,fc,fb,Q0) tries to % approximate the integral of f(x) from a to b to within a % relative error of tol. F is a string containing the name % of f. The remaining arguments are generated by quad or by % the recursion. lev is the recursion level. % fa = f(a). fc = f((a+b)/2). fb = f(b). % Q0 is an approximate value of the integral. % See also QUAD and QUAD8. % C.B. Moler, 3-22-87. LEVMAX = 100; if lev > LEVMAX % Print warning the first time the recursion level is exceeded. if ~recur_lev_excess % warning('Recursion level limit reached in quad. Singularity likely.') recur_lev_excess = 1; else recur_lev_excess = recur_lev_excess + 1; end Q = Q0; cnt = 0; c = (a + b)/2; if trace yc = feval(FunFcn,c,varargin{:}); plot(c,yc,'.'); if ~isreal(yc) plot(c,imag(yc),'+'); end drawnow end else % Evaluate function at midpoints of left and right half intervals. h = b - a; c = (a + b)/2; x = [a+h/4, b-h/4]; f = feval(FunFcn,x,varargin{:}); if trace plot(x,f,'.'); if ~isreal(f) plot(x,imag(f),'+') end drawnow end cnt = 2; % Simpson's rule for half intervals. Q1 = h*(fa + 4*f(1) + fc)/12; Q2 = h*(fc + 4*f(2) + fb)/12; Q = Q1 + Q2; % Recursively refine approximations. if abs(Q - Q0) > tol(1)*abs(Q) + tol(2); [Q1,cnt1,recur_lev_excess] = quadstp(FunFcn,a,c,tol/2,lev+1,fa,f(1),fc,Q1,trace,recur_lev_excess,varargin{:}); [Q2,cnt2,recur_lev_excess] = quadstp(FunFcn,c,b,tol/2,lev+1,fc,f(2),fb,Q2,trace,recur_lev_excess,varargin{:}); Q = Q1 + Q2; cnt = cnt + cnt1 + cnt2; end end