invsi.m

很多matlab的源代码

function [b,k]=invsi(y,a,b,tol)
% INVSI Inverse of sine integral with argument at.
%
%	[C,K] = INVSI(X,A,B,TOL) The inverse of the sine integral 
%	If X = SI(Y,A) or the integral of sin(At)/(At) over [0, Y]
%	Then C = INVSI(X,A) should equal Y to within approximation error.
%	A = argument in sin(At)/(At) [Default: A = 1].
%       NOTE: Value C returned is the SMALLEST possible choice over (0,pi/A).
%
%	B = initial guess. The value C may change with the starting guess B
%	NOTE: B must be a scalar OR the same size as X [Default: B=0].
%	TOL = desired tolerance [Default: TOL = 1e-7]
%       K returns the number of iterations required.
%
%       INVSI (with no input arguments) invokes the following example:
%
%	% Let Y0 = [0.6 2.4], Ys = si(Y0). Then X0=invsi(Ys) should equal Y0 
%	  >>Y0 = [0.6 2.4], Ys = si(Y0), X0 = invsi(Ys) 
%       % BUT if Y1 = 4, Ys1 = si(Y1), Then X1 = invsi(Ys1) does not equal 4 
%         >>Y1=4,Ys1=si(Y1),X1 = invsi(Ys1) 
%       % X1 = 2.4206 since si(4) = si(2.4206) = 1.7582
%	% If we use an initial guess of 5, X2=invsi(Ys1,1.2) does equal 4
%         >>invsi(Ys1,1,5)  %Note: a = 1 and b = 5


% 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 invsi,disp('Strike a key to see results of the example')
pause,Y0 = [0.6 2.4], Ys = si(Y0), X0 = invsi(Ys),disp('Hit a key to continue')
pause,Y1=4,Ys1=si(Y1),X1 = invsi(Ys1),b = invsi(Ys1,1,5);return,end

if nargin<4,tol=1e-7;end,
if nargin<3,b=0;end,
if nargin<2,a=1;end,
z1=si(pi);
if prod(size(b))~=1 & size(y)~=size(b),error('Inconsistent dimensions'),end
z=abs(a*y);x=z.*(z<z1);[m,n]=size(x);
bp=ones(m,n);b=b.*bp;
kmax=100; %Change as required
k=0;
while any(any(abs(b-bp)>tol)) & k<=kmax
k=k+1;bp=b;f=(si(bp)-x)/2;
fd=sin(bp+eps)./(bp+eps);b=bp-f./fd;
end
if k>kmax,disp(['No convergence in ' num2str(kmax) ' iterations']),end
b=b.*sign(y)/abs(a);
if any(any(z>z1)),i=find(z>z1);[m,n]=size(i);b(i)=nan*ones(m,n);end