Необходимые программы-функции
1. Обратное преобразование Фурье
function x = ifft2(f, mrows, ncols)
%IFFT2 Two-dimensional inverse discrete Fourier transform.
% IFFT2(F) returns the two-dimensional inverse Fourier transform
% of matrix F. If F is a vector, the result will have the same
% orientation.
%
% IFFT2(F,MROWS,NCOLS) pads matrix F with zeros to size MROWS-by-NCOLS
% before transforming.
%
% See also FFT2, IFFT, FFTSHIFT.
% J.N. Little 12-18-85
% Revised 4-15-87 JNL
% Revised 5-3-90 CRD
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 5.14 $ $Date: 2000/08/15 16:49:31 $
if (~isa(f, 'double'))
f = double(f);
end
if ndims(f)==2
if nargin==1
x = ifftn(f);
else
x = ifftn(f,[mrows ncols]);
end
else
f = conj(f);
mn = size(f,1)*size(f,2);
if nargin==1
x = conj(fft(fft(f,[],2),[],1))/mn;
else
x = conj(fft(fft(f,ncols,2),mrows,1))/(mrows*ncols);
end
end
2.N-мерное обратное преобразование Фурье
%IFFTN N-dimensional inverse discrete Fourier transform.
% IFFTN(F) returns the N-dimensional inverse discrete Fourier
% transform of the N-D array F. If F is a vector, the result
% will have the same orientation.
%
% IFFTN(F,SIZ) pads F so that its size vector is SIZ before
% performing the transform. If any element of SIZ is smaller
% than the corresponding dimension of F, then F will be cropped
% in that dimension.
%
% See also FFT, IFFT, FFT2, IFFT2, FFTN, FFTSHIFT.
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 1.12 $ $Date: 2000/06/01 02:10:46 $
3. Прямое преобразование Фурье
%FFT Discrete Fourier transform.
% FFT(X) is the discrete Fourier transform (DFT) of vector X. For
% matrices, the FFT operation is applied to each column. For N-D
% arrays, the FFT operation operates on the first non-singleton
% dimension.
%
% FFT(X,N) is the N-point FFT, padded with zeros if X has less
% than N points and truncated if it has more.
%
% FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
% dimension DIM.
%
% For length N input vector x, the DFT is a length N vector X,
% with elements
% N
% X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
% n=1
% The inverse DFT (computed by IFFT) is given by
% N
% x(n) = (1/N) sum X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
% k=1
%
% The relationship between the DFT and the Fourier coefficients a and b in
% N/2
% x(n) = a0 + sum a(k)*cos(2*pi*k*t(n)/(N*dt))+b(k)*sin(2*pi*k*t(n)/(N*dt))
% k=1
% is
% a0 = X(1)/N, a(k) = 2*real(X(k+1))/N, b(k) = -2*imag(X(k+1))/N,
% where x is a length N discrete signal sampled at times t with spacing dt.
%
% See also IFFT, FFT2, IFFT2, FFTSHIFT.
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 5.12 $ $Date: 2000/06/14 21:05:05 $
% Built-in function.
4. Объединение комплексных векторов
%CONJ Complex conjugate.
% CONJ(X) is the complex conjugate of X.
% For a complex X, CONJ(X) = REAL(X) - i*IMAG(X).
%
% See also REAL, IMAG, I, J.
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 5.6 $ $Date: 2000/06/01 00:40:15 $
% Built-in function.
5. ФУНКЦИИ logspace и freqspace
Generate logarithmically spaced vectors
Syntax
y = logspace(a,b)
y = logspace(a,b,n)
y = logspace(a,pi)
Description
The logspace function generates logarithmically spaced vectors. Especially useful for
creating frequency vectors, it is a logarithmic equivalent of linspace and the ":" or
colon operator.
y = logspace(a,b) generates a row vector y of 50 logarithmically spaced points
between decades 10^a and 10^b.
y = logspace(a,b,n) generates n points between decades 10^a and 10^b.
y = logspace(a,pi) generates the points between 10^a and pi, which is useful for
digital signal processing where frequencies over this interval go around the unit circle.
Remarks
All the arguments to logspace must be scalars.
freqspace
Frequency spacing for frequency response.
Syntax
f = freqspace(n)
f = freqspace(n,'whole')
[f1,f2] = freqspace(n)
[f1,f2] = freqspace([m n])
[x1,y1] = freqspace(n,'meshgrid')
[x1,y1] = freqspace([m n],'meshgrid')
Description
freqspace returns the implied frequency range for equally spaced frequency
responses. This is useful when creating frequency vectors for use with freqz.
f = freqspace(n) returns the frequency vector f assuming n evenly spaced points
around the unit circle. For n even or odd, f is (0:2/n:1). For n even, freqspace
returns (n + 2)/2 points. For N odd, it returns (n + 1)/2 points.
f = freqspace(n,'whole') returns n evenly spaced points around the whole unit
circle. In this case, f is 0:2/n:2*(n-1)/n.
[f1,f2] = freqspace(n) returns the two-dimensional frequency vectors f1 and f2
for an n-by-n matrix. For n odd, both f1 and f2 are [-1 + 1/n:2/n:1-1/n]. For n
even, both f1 and f2 are [-1:2/n:1-2/n].
[f1,f2] = freqspace([m n]) returns the two-dimensional frequency vectors f1
and f2 for an m-by-n matrix.
[x1,y1] = freqspace(n,'meshgrid') and
[x1,y1] = freqspace([m n],'meshgrid') are equivalent to
[f1,f2] = freqspace(...);
[x1,y1] = meshgrid(f1,f2);
See the MATLAB documentation for details on the meshgrid function.
6. Преобразование Гильберта
hilbert
Compute the discrete-time analytic signal using the Hilbert transform.
Syntax
x = hilbert(xr)
x = hilbert(xr,n)
Description
x = hilbert(xr) returns a complex helical sequence, sometimes called the analytic
signal, from a real data sequence. The analytic signal x = xr + i*xi has a real part,
xr, which is the original data, and an imaginary part, xi, which contains the Hilbert
transform. The imaginary part is a version of the original real sequence with a 90°
phase shift. Sines are therefore transformed to cosines and vice versa. The Hilbert
transformed series has the same amplitude and frequency content as the original real
data and includes phase information that depends on the phase of the original data.
If xr is a matrix, x = hilbert(xr) operates columnwise on the matrix, finding the
Hilbert transform of each column.
x = hilbert(xr,n) uses an n point FFT to compute the Hilbert transform. The input
data xr is zero-padded or truncated to length n, as appropriate.
The Hilbert transform is useful in calculating instantaneous attributes of a time series,
especially the amplitude and frequency. The instantaneous amplitude is the amplitude
of the complex Hilbert transform; the instantaneous frequency is the time rate of change
of the instantaneous phase angle. For a pure sinusoid, the instantaneous amplitude
and frequency are constant. The instantaneous phase, however, is a sawtooth,
reflecting the way in which the local phase angle varies linearly over a single cycle. For
mixtures of sinusoids, the attributes are short term, or local, averages spanning no
more than two or three points.
Reference [1] describes the Kolmogorov method for minimum phase reconstruction,
which involves taking the Hilbert transform of the logarithm of the spectral density of a
time series. The toolbox function rceps performs this reconstruction.
For a discrete-time analytic signal x, the last half of fft(x) is zero, and the first (DC)
and center (Nyquist) elements of fft(x) are purely real.
Algorithm
The analytic signal for a sequence x has a one-sided Fourier transform, that is,
negative frequencies are 0. To approximate the analytic signal, hilbert calculates the
FFT of the input sequence, replaces those FFT coefficients that correspond to negative
frequencies with zeros, and calculates the inverse FFT of the result.
In detail, hilbert uses a four-step algorithm:
It calculates the FFT of the input sequence, storing the result in a vector x.
It creates a vector h whose elements h(i) have the values:
1 for i = 1, (n/2)+1
2 for i = 2, 3, ... , (n/2)
0 for i = (n/2)+2, ... , n
It calculates the element-wise product of x and h.
It calculates the inverse FFT of the sequence obtained in step 3 and returns the
first n elements of the result.
If the input data xr is a matrix, hilbert operates in a similar manner, extending each
step above to handle the matrix case.
Example
xr = [1 2 3 4];
x = hilbert(xr)
x =
1.0000+1.0000i 2.0000-1.0000i 3.0000-1.0000i 4.0000+1.0000i
You can see that the imaginary part, imag(x) = [1 -1 -1 1], is the Hilbert transform
of xr, and the real part, real(x) = [1 2 3 4], is simply xr itself. Note that the last
half of fft(x) = [10 -4+4i -2 0] is zero (in this example, the last half is just the last
element), and that the DC and Nyquist elements of fft(x), 10 and -2 respectively, are
purely real.
See Also
rceps
Real cepstrum and minimum phase reconstruction.
References
[1] Claerbout, J.F., Fundamentals of Geophysical Data Processing, McGraw-Hill, 1976,
pp. 59-62.
[2] Marple, S.L., "Computing the discrete-time analytic signal via FFT," IEEE
Transactions on Signal Processing, Vol. 47, No. 9 (September 1999), pp. 2600-2603.
[3] Oppenheim, A.V., and R.W. Schafer, Discrete-Time Signal Processing, 2nd ed.,
Prentice-Hall, 1998.
7. Ряд Тейлора
function E = expm2(A)
%EXPM2 Matrix exponential via Taylor series.
% E = expm2(A) illustrates the classic definition for the
% matrix exponential. As a practical numerical method,
% this is often slow and inaccurate.
%
% See also EXPM, EXPM1, EXPM3.
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 5.6 $ $Date: 2000/06/01 02:04:14 $
% Taylor series for exp(A)
E = zeros(size(A));
F = eye(size(A));
k = 1;
while norm(E+F-E,1) > 0
E = E + F;
F = A*F/k;
k = k+1;
end
%EIG Eigenvalues and eigenvectors.
% E = EIG(X) is a vector containing the eigenvalues of a square
% matrix X.
%
% [V,D] = EIG(X) produces a diagonal matrix D of eigenvalues and a
% full matrix V whose columns are the corresponding eigenvectors so
% that X*V = V*D.
%
% [V,D] = EIG(X,'nobalance') performs the computation with balancing
% disabled, which sometimes gives more accurate results for certain
% problems with unusual scaling.
%
% E = EIG(A,B) is a vector containing the generalized eigenvalues
% of square matrices A and B.
%
% [V,D] = EIG(A,B) produces a diagonal matrix D of generalized
% eigenvalues and a full matrix V whose columns are the
% corresponding eigenvectors so that A*V = B*V*D.
%
% EIG(A,B,'chol') is the same as EIG(A,B) for symmetric A and symmetric
% positive definite B. It computes the generalized eigenvalues of A and B
% using the Cholesky factorization of B.
% EIG(A,B,'qz') ignores the symmetry of A and B and uses the QZ algorithm.
% In general, the two algortihms return the same result, however using the
% QZ algorithm may be more stable for certain problems.
% The flag is ignored when A and B are not symmetric.
%
% See also CONDEIG, EIGS.
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 5.9 $ $Date: 2000/06/01 02:04:14 $
% Built-in function.
8. Алгоритм-функция определения нулей и полюсов
function [z,p,k] = tf2zp(num,den)
%TF2ZP Transfer function to zero-pole conversion.
% [Z,P,K] = TF2ZP(NUM,DEN) finds the zeros, poles, and gains:
%
% (s-z1)(s-z2)...(s-zn)
% H(s) = K ---------------------
% (s-p1)(s-p2)...(s-pn)
%
% from a SIMO transfer function in polynomial form:
%
% NUM(s)
% H(s) = --------
% DEN(s)
%
% Vector DEN specifies the coefficients of the denominator in
% descending powers of s. Matrix NUM indicates the numerator
% coefficients with as many rows as there are outputs. The zero
% locations are returned in the columns of matrix Z, with as many
% columns as there are rows in NUM. The pole locations are returned
% in column vector P, and the gains for each numerator transfer
% function in vector K.
% For discrete-time transfer functions, it is highly recommended to
% make the length of the numerator and denominator equal to ensure
% correct results. You can do this using the function EQTFLENGTH in
% the Signal Processing Toolbox. However, this function only handles
% single-input single-output systems.
% See also ZP2TF.
% Clay M. Thompson 11-6-90
% Revised ACWG 1-17-90
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 1.24 $ $Date: 2000/01/31 04:12:17 $
[num,den] = tfchk(num,den);
% Normalize transfer function
if length(den)
coef = den(1);
else
coef = 1;
end
if abs(coef)<eps,
error('Denominator must have non-zero leading coefficient.');
end
den = den./coef;
num = num./coef;
% Remove leading columns of zeros from numerator
if length(num)
while(all(num(:,1)==0) & length(num) > 1)
num(:,1) = [];
end
end
[ny,np] = size(num);
% Poles
p = roots(den);
% Zeros and Gain
k = zeros(ny,1);
linf = inf;
z = linf(ones(np-1,1),ones(ny,1));
for i=1:ny
zz = roots(num(i,:));
if length(zz), z(1:length(zz), i) = zz; end
ndx = find(num(i,:)~=0);
if length(ndx), k(i,1) = num(i,ndx(1)); end
end
9. Функция- интерполянт рядом Фурье
function y = interpft(x,ny,dim);
%INTERPFT 1-D interpolation using FFT method.
% Y = INTERPFT(X,N) returns a vector Y of length N obtained
% by interpolation in the Fourier transform of X.
% If X is a matrix, interpolation is done on each column.
% If X is an array, interpolation is performed along the first
% non-singleton dimension.
% INTERPFT(X,N,DIM) performs the interpolation along the
% dimension DIM.
% Assume x(t) is a periodic function of t with period p, sampled
% at equally spaced points, X(i) = x(T(i)) where T(i) = (i-1)*p/M,
% i = 1:M, M = length(X). Then y(t) is another periodic function
% with the same period and Y(j) = y(T(j)) where T(j) = (j-1)*p/N,
% j = 1:N, N = length(Y). If N is an integer multiple of M,
% then Y(1:N/M:N) = X.
% See also INTERP1.
% Charles R. Denham, MathWorks, 1988.
% Robert Piche, Tampere University of Technology, 10/93.
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 5.13 $ $Date: 2000/06/02 17:51:32 $
error(nargchk(2,3,nargin));
if nargin==2,
[x,nshifts] = shiftdim(x);
if prod(size(x))==1, nshifts = 1; end % Return a row for a scalar
elseif nargin==3,
perm = [dim:max(length(size(x)),dim) 1:dim-1];
x = permute(x,perm);
end
siz = size(x);
[m,n] = size(x);
if prod(size(ny))~=1, error('N must be a scalar.'); end
% If necessary, increase ny by an integer multiple to make ny > m.
if ny > m
incr = 1;
else
if ny==0, y=[]; return, end
incr = floor(m/ny) + 1;
ny = incr*ny;
end
a = fft(x,[],1);
nyqst = ceil((m+1)/2);
b = [a(1:nyqst,:) ; zeros(ny-m,n) ; a(nyqst+1:m,:)];
if rem(m,2) == 0
b(nyqst,:) = b(nyqst,:)/2;
b(nyqst+ny-m,:) = b(nyqst,:);
end
y = ifft(b,[],1);
if isreal(x), y = real(y); end
y = y * ny / m;
y = y(1:incr:ny,:); % Skip over extra points when oldny <= m.
if nargin==2,
y = reshape(y,[ones(1,nshifts) size(y,1) siz(2:end)]);
elseif nargin==3,
y = ipermute(y,perm);
end
10. N-мерный интерполянт
function vi = interpn(varargin)
%INTERPN N-D interpolation (table lookup).
% VI = INTERPN(X1,X2,X3,...,V,Y1,Y2,Y3,...) interpolates to find VI,
% the values of the underlying N-D function V at the points in
% arrays Y1,Y2,Y3,etc. For an N-D V, INTERPN should be called with
% 2*N+1 arguments. Arrays X1,X2,X3,etc. specify the points at which
% the data V is given. Out of range values are returned as NaN's.
% Y1,Y2,Y3,etc. must be arrays of the same size or vectors. Vector
% arguments that are not the same size, and have mixed orientations
% (i.e. with both row and column vectors) are passed through NDGRID to
% create the Y1,Y2,Y3,etc. arrays. INTERPN works for all N-D arrays
% with 2 or more dimensions.
% VI = INTERPN(V,Y1,Y2,Y3,...) assumes X1=1:SIZE(V,1),X2=1:SIZE(V,2),etc.
% VI = INTERPN(V,NTIMES) expands V by interleaving interpolates between
% every element, working recursively for NTIMES.
% VI = INTERPN(V) is the same as INTERPN(V,1).
% VI = INTERPN(...,'method') specifies alternate methods. The default
% is linear interpolation. Available methods are:
% 'linear' - linear interpolation
% 'cubic' - cubic interpolation
% 'nearest' - nearest neighbor interpolation
% 'spline' - spline interpolation
% INTERPN requires that X1,X2,X3,etc. be monotonic and plaid (as if
% they were created using NDGRID). X1,X2,X3,etc. can be non-uniformly
% spaced. For faster interpolation when X1,X2,Y3,etc. are known to be
% equally spaced and monotonic, use the methods '*linear', '*cubic', or
% '*nearest'.
% See also INTERP1, INTERP2, INTERP3, NDGRID.
% Clay M. Thompson 8-1-94
% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 1.30 $ $Date: 2000/06/08 19:57:22 $
bypass = 0;
uniform = 1;
if isstr(varargin{end}),
narg = nargin-1;
method = [varargin{end} ' ']; % Protect against short string.
if method(1)=='*', % Direct call bypass.
if method(2)=='l', % linear interpolation.
vi = linear(varargin{1:end-1});
return
elseif method(2)=='c', % cubic interpolation
vi = cubic(varargin{1:end-1});
return
elseif method(2)=='n', % Nearest neighbor interpolation
vi = nearest(varargin{1:end-1});
return
elseif method(2)=='s', % spline interpolation
method = 'spline'; bypass = 1;
else
error([deblank(method),' is an invalid method.']);
end
elseif method(1)=='s', % Spline interpolation
method = 'spline'; bypass = 1;
end
else
narg = nargin;
method = 'linear';
end
if narg==0,
error('Not enough input arguments.');
elseif narg<=2. % interp(v) or interpn(v,n), Expand V ntimes
if narg==1, ntimes = 1; else ntimes = floor(varargin{2}); end
siz = size(varargin{1});
x = cell(size(siz)); y = cell(size(siz));
for i=1:length(siz),
x{i} = 1:siz(i);
y{i} = 1:1/(2^ntimes):siz(i);
end
[msg,x,v,y] = xnchk(x,varargin{1},y);
elseif all(size(varargin{1}) >= 2) & (narg == ndims(varargin{1})+1),
% If every dimension of V has size >=2 and nargin = ndims + 1,
% then handle interpn(v,y1,y2,y3,...).
v = varargin{1};
if ndims(v)~=narg-1, error('Wrong number of input arguments.'); end
siz = size(v);
y = varargin(2:narg);
x = cell(size(siz));
for i=1:length(siz), x{i} = 1:siz(i); end
[msg,x,v,y] = xnchk(x,v,y);
elseif (rem(narg,2)==1) & all(size(varargin{floor((narg+1)/2)}) >= 2) & ...
(narg == 2*ndims(varargin{(narg+1)/2})+1),
% If the number of input arguments is odd and every dimension of V
% has size >= 2 and nargin = 2*ndims+1,
% then handle interpn(x1,x2,x3,...,v,y1,y2,y3,...)
v = varargin{(narg+1)/2};
siz = size(v);
x = varargin(1:ndims(v));
y = varargin(ndims(v)+2:narg);
[msg,x,v,y] = xnchk(x,v,y);
else
error(...
'Wrong number of input arguments or some dimension of V is less than 2.');
end
if ~isempty(msg), error(msg); end
if ~bypass,
% Create xx cell array containing the vectors from each Xi array
xx = cell(size(x));
for i=1:length(x),
ind(1:length(x)) = {1};
ind{i} = ':';
xx{i} = x{i}(ind{:});
xx{i} = xx{i}(:); % Make sure its a column
end
% Check for non-equally spaced data. If so, map(x1,x2,x3,...) and
% (y1,y2,y3,...) to array index coordinate system.
for i=1:prod(size(xx));
dx{i} = diff(xx{i});
xxmax(i) = max(abs(xx{i}));
if prod(size(dx{i}))==1,
maxabs(i) = 0;
else
maxabs(i) = max(abs(diff(dx{i})));
end
end
if any(maxabs > eps*xxmax), % If data is not equally spaced,
for i=1:length(x),
% Flip orientation of data so that x{i} is increasing
if any(dx{i} < 0),
for j=1:length(x), x{j} = flipdim(x{j},i); end
for j=1:length(y), y{j} = flipdim(y{j},i); end
v = flipdim(v,i);
xx{i} = flipud(xx{i});
dx{i} = -flipud(dx{i});
end
end
for i=1:length(dx),
if any(dx{i}<=0),
error('X1,X2,X3, etc. must be monotonic vectors or arrays produced by NDGRID.');
end
end
% Bypass mapping code for cubic
if method(1)~='c',
% Map values in y to values in si via linear interpolation (one for
% each input argument.
for k=1:length(y),
zz = xx{k}; zi = y{k}; dz = dx{k}; si = [];
% Determine the nearest location of zi in zz
[zzi,j] = sort(zi(:));
[dum,i] = sort([zz;zzi]);
si(i) = (1:length(i));
si = (si(length(zz)+1:end)-(1:length(zzi)))';
si(j) = si;
% Map values in zi to index offset (si) tia linear interpolation
si(si<1) = 1;
si(si>length(zz)-1) = length(zz)-1;
si = si + (zi(:)-zz(si))./(zz(si+1)-zz(si));
y{k}(:) = si; % Replace value in y
end
% Change x1,x2,x3,... to be index cordinates.
for i=1:length(x),
x{i} = 1:siz(i);
end
[x{1:end}] = ndgrid(x{:});
else
uniform = 0;
end
end
end
% Now do the interpolation based on method.
method = [lower(method),' ']; % Protect against short string
if method(1)=='l', % linear interpolation.
vi = linear(x{:},v,y{:});
elseif method(1)=='c', % cubic interpolation
if uniform
vi = cubic(x{:},v,y{:});
else
d = zeros(size(y{1}));
for i=1:length(y)
d = d | y{i} < min(x{i}(:)) | y{i} > max(x{i}(:));
end
d = find(d);
vi = splinen(x,v,y);
vi(d) = NaN;
end
elseif method(1)=='n', % Nearest neighbor interpolation
vi = nearest(x{:},v,y{:});
elseif method(1)=='s', % Spline interpolation
vi = splinen(x,v,y);
else
error([deblank(method),' is an invalid method.']);
end
%------------------------------------------------------
function F = linear(varargin)
%LINEAR N-D linear data interpolation.
% VI = LINEAR(X1,X2,X3,...,V,Y1,Y2,Y3,...) interpolates to find VI,
% the values of the underlying N-D function V, at the points in arrays
% Y1,Y2,Y3, etc. via linear interpolation. For an N-D V, LINEAR
% should be called with 2*N+1 arguments. Arrays X1,X2,X3,... specify
% the points at which the data V is given. X1,X2,X3, etc. can also be
% vectors specifying the abscissae for the matrix V as for NDGRID. In
% both cases, X1,X2,X3, etc. must be equally spaced and monotonic.
% Out of range values are returned as NaN.
%
% VI = LINEAR(V,Y1,Y2,Y3,...) assumes X1=1:SIZE(V,1), X2=1:SIZE(V,2),
% X3=1:SIZE(V,3), etc.
% VI = LINEAR(V,NTIMES) expands V by interleaving interpolates between
% every element, working recursively for NTIMES. LINEAR(V) is the
% same as LINEAR(V,1).
%
% See also INTERPN.
% Clay M. Thompson 8-2-94
if nargin==0,
error('Not enough input arguments.');
elseif nargin<=2, % linear(v) or linear(v,n), Expand V n times
if nargin==1, ntimes = 1; else ntimes = floor(varargin{2}); end
v = varargin{1};
siz = size(v);
s = cell(size(siz));
for i=1:length(s),
s{i} = 1:1/(2^ntimes):siz(i);
end
[s{1:end}] = ndgrid(s{:});
elseif nargin>2 & (rem(nargin,2)==0 | nargin == ndims(varargin{1})+1),
% linear(v,y1,y2,y3,...)
v = varargin{1};
if ndims(v)~=nargin-1,
error('Wrong number of input arguments.');
end
siz = size(v);
s = varargin(2:nargin);
elseif nargin>2 & rem(nargin,2)==1 & ...
(nargin == 2*ndims(varargin{(nargin+1)/2})+1),
% linear(x1,x2,x3,...,v,y1,y2,y3,...)
v = varargin{(nargin+1)/2};
if nargin~=2*ndims(v)+1, error('Wrong number of input arguments.'); end
siz = size(v);
x = varargin(1:ndims(v));
y = varargin(ndims(v)+2:nargin);
[msg,x,v,y] = xnchk(x,v,y);
for i=length(x):-1:1,
xsiz{i} = size(x{i});
xprodsiz(i) = prod(size(x{i}));
end
if ~isequal(xprodsiz,siz) & ~isequal(siz,xsiz{:}),
error('The lengths of the X1,X2,X3,... vectors must match V.');
end
s = cell(size(y));
for i=1:length(s),
s{i} = 1 + (y{i}-x{i}(1))/(x{i}(xprodsiz(i))-x{i}(1))*(siz(i)-1);
end
end
if any(siz<2*ones(size(siz))),
error('V must be at least 2-by-2-by-2-by-...');
end
for i=length(s):-1:1, ssiz{i} = size(s{i}); end
if ~isequal(ssiz{:}),
error('Y1,Y2,Y3, etc. must be the same size.');
end
% Check for out of range values of s and set to 1
sout = cell(size(s));
for i=1:length(s),
sout{i} = find((s{i}<1) | (s{i}>siz(i)));
if ~isempty(sout{i}), s{i}(sout{i}) = ones(size(sout{i})); end
end
% Matrix element indexing
offset = cumprod([1 siz(1:end-1)]);
ndx = 1;
for i=1:length(s),
ndx = ndx + offset(i)*floor(s{i}-1);
end
% Compute intepolation parameters, check for boundary value.
for i=1:length(s),
if isempty(s{i}), d = s{i}; else d = find(s{i}==siz(i)); end
s{i} = s{i}-floor(s{i});
if ~isempty(d), s{i}(d) = s{i}(d)+1; ndx(d) = ndx(d)-offset(i); end
end
d = []; % Reclaim memory
% Create index arrays, iw.
iw = cell(size(s));
[iw{1:end}] = ndgrid(0:1);
% Reshape each iw{i} to a column and then arrange into a matrix
iwcol = [prod(size(iw{1})) 1];
for i=1:prod(size(iw)), iw{i} = reshape(iw{i},iwcol); end
iw = cat(2,iw{:}); % Arrange columns into a matrix
% Do the linear interpolation: f = v1*(1-s) + v2*s along each direction
F = 0;
for i=1:size(iw,1),
vv = v(ndx + offset*iw(i,:)');
for j=1:size(iw,2),
switch iw(i,j)
case 0 % Interpolation function (1-s)
vv = vv.*(1-s{j});
case 1 % Interpolation function (s)
vv = vv.*s{j};
end
end
F = F + vv;
end
% Now set out of range values to NaN.
for i=1:length(sout),
if ~isempty(sout{i}), F(sout{i}) = NaN; end
end
%------------------------------------------------------
function F = cubic(varargin)
%CUBIC Cubic data interpolation.
% CUBIC(...) is the same as LINEAR(....) except that it uses
% cubic interpolation.
%
% See also INTERPN.
% Clay M. Thompson 8-1-94
% Based on "Cubic Convolution Interpolation for Digital Image
% Processing", Robert G. Keys, IEEE Trans. on Acoustics, Speech, and
% Signal Processing, Vol. 29, No. 6, Dec. 1981, pp. 1153-1160.
if nargin==0,
error('Not enough input arguments.');
elseif nargin<=2, % cubic(v) or cubic(v,n), Expand V n times
if nargin==1, ntimes = 1; else ntimes = floor(varargin{2}); end
v = varargin{1};
siz = size(v);
s = cell(size(siz));
for i=1:length(s),
s{i} = 1:1/(2^ntimes):siz(i);
end
[s{1:end}] = ndgrid(s{:});
elseif nargin>2 & (rem(nargin,2)==0 | nargin == ndims(varargin{1})+1),
% cubic(v,y1,y2,y3,...)
v = varargin{1};
if ndims(v)~=nargin-1,
error('Wrong number of input arguments.');
end
siz = size(v);
s = varargin(2:nargin);
elseif nargin>2 & rem(nargin,2)==1 & ...
(nargin == 2*ndims(varargin{(nargin+1)/2})+1),
% cubic(x1,x2,x3,...,v,y1,y2,y3,...)
v = varargin{(nargin+1)/2};
if nargin~=2*ndims(v)+1, error('Wrong number of input arguments.'); end
siz = size(v);
x = varargin(1:ndims(v));
y = varargin(ndims(v)+2:nargin);
[msg,x,v,y] = xnchk(x,v,y);
for i=length(x):-1:1,
xsiz{i} = size(x{i});
xprodsiz(i) = prod(size(x{i}));
end
if ~isequal(xprodsiz,siz) & ~isequal(siz,xsiz{:}),
error('The lengths of the X1,X2,X3,... vectors must match V.');
end
s = cell(size(y));
for i=1:length(s),
s{i} = 1 + (y{i}-x{i}(1))/(x{i}(xprodsiz(i))-x{i}(1))*(siz(i)-1);
end
end
if any(siz<3*ones(size(siz))),
error('V must be at least 3-by-3-by-3-by-...');
end
for i=length(s):-1:1, ssiz{i} = size(s{i}); end
if ~isequal(ssiz{:}),
error('Y1,Y2,Y3, etc. must be the same size.');
end
% Check for out of range values of s and set to 1
sout = cell(size(s));
for i=1:length(s),
sout{i} = find((s{i}<1) | (s{i}>siz(i)));
if ~isempty(sout{i}), s{i}(sout{i}) = ones(size(sout{i})); end
end
% Matrix element indexing
offset = cumprod([1 siz(1:end-1)+2]);
ndx = 1;
for i=1:length(s),
ndx = ndx + offset(i)*floor(s{i}-1);
end
% Compute intepolation parameters, check for boundary value.
for i=1:length(s),
if isempty(s{i}), d = s{i}; else d = find(s{i}==siz(i)); end
s{i} = s{i}-floor(s{i});
if ~isempty(d), s{i}(d) = s{i}(d)+1; ndx(d) = ndx(d)-offset(i); end
end
d = []; % Reclaim memory
%
% Expand v so interpolation is valid at the boundaries.
%
vv = zeros(siz+2);
for i=length(siz):-1:1,
ind{i} = 2:siz(i)+1;
end
vv(ind{:}) = v;
% Insert values on the boundary that allow the use of the same
% interpolation kernel everywhere. These values are computed
% computed from the 3 strips nearest the boundary:
% v(boundary,:) = 3*v(1,:) - 3*v(2,:) + v(3,:);
% Use comma separated list equivalence to index into vv. That is,
% vv(edges{1,:}) is the same as vv(edges{1,1},edges{1,2},edges{1,3},...).
for i=length(siz):-1:1,
ind{i} = 1:siz(i)+2; % vv(ind{:}) is the same as vv(:,:,...,:).
end
for i=1:length(s),
edges = ind(ones(4,1),:); % Reinitialize edges
% Set virtual right edge
edges(:,i) = {1;2;3;4};
vv(edges{1,:}) = 3*vv(edges{2,:}) - 3*vv(edges{3,:}) + vv(edges{4,:});
% Set virtual left edge
edges(:,i) = num2cell(siz(i)+[2 1 0 -1]');
vv(edges{1,:}) = 3*vv(edges{2,:}) - 3*vv(edges{3,:}) + vv(edges{4,:});
end
siz = siz + 2;
% Create index arrays, iw.
iw = cell(size(s));
[iw{1:end}] = ndgrid(0:3);
% Reshape each iw{i} to a column and then arrange into a matrix
iwcol = [prod(size(iw{1})) 1];
for i=1:prod(size(iw)), iw{i} = reshape(iw{i},iwcol); end
iw = cat(2,iw{:});
% Do the cubic interpolation:
% f = v1*((2-s)*s-1)*s + v2*((3*s-5)*s^2+2) + ...
% v3*(((4-3*s)*s+1)*s) + v4*((s-1)*s^2) along each direction
F = 0;
for i=1:size(iw,1),
vvv = vv(ndx + offset*iw(i,:)');
for j=1:size(iw,2),
switch iw(i,j)
case 0 % Interpolation function ((2-s)*s-1)*s
vvv = vvv.*( ((2-s{j}).*s{j}-1).*s{j} );
case 1 % Interpolation function (3*s-5)*s*s+2
vvv = vvv.*( (3*s{j}-5).*s{j}.*s{j}+2 );
case 2 % Interpolation function ((4-3*s)*s+1)*s
vvv = vvv.*( ((4-3*s{j}).*s{j}+1).*s{j} );
case 3 % Interpolation function (s-1)*s*s
vvv = vvv.*( (s{j}-1).*s{j}.*s{j} );
end
end
F = F + vvv;
end
F(:) = F/pow2(length(s));
% Now set out of range values to NaN.
for i=1:length(sout),
if ~isempty(sout{i}), F(sout{i}) = NaN; end
end
%------------------------------------------------------
function F = nearest(varargin)
%NEAREST 3-D Nearest neighbor interpolation.
% NEAREST(...) is the same as LINEAR(...) except that it uses
% nearest neighbor interpolation.
%
% See also INTERPN.
% Clay M. Thompson 1-31-94
if nargin==0,
error('Not enough input arguments.');
elseif nargin<=2, % nearest(v) or nearest(v,n), Expand V n times
if nargin==1, ntimes = 1; else ntimes = floor(varargin{2}); end
v = varargin{1};
siz = size(v);
s = cell(size(siz));
for i=1:length(s),
s{i} = 1:1/(2^ntimes):siz(i);
end
[s{1:end}] = ndgrid(s{:});
elseif nargin>2 & (rem(nargin,2)==0 | nargin == ndims(varargin{1})+1),
% nearest(v,y1,y2,y3,...)
v = varargin{1};
if ndims(v)~=nargin-1,
error('Wrong number of input arguments.');
end
siz = size(v);
s = varargin(2:nargin);
elseif nargin>2 & rem(nargin,2)==1 & ...
(nargin == 2*ndims(varargin{(nargin+1)/2})+1),
% nearest(x1,x2,x3,...,v,y1,y2,y3,...)
v = varargin{(nargin+1)/2};
if nargin~=2*ndims(v)+1, error('Wrong number of input arguments.'); end
siz = size(v);
x = varargin(1:ndims(v));
y = varargin(ndims(v)+2:nargin);
[msg,x,v,y] = xnchk(x,v,y);
for i=length(x):-1:1,
xsiz{i} = size(x{i});
xprodsiz(i) = prod(size(x{i}));
end
if ~isequal(xprodsiz,siz) & ~isequal(siz,xsiz{:}),
error('The lengths of the X1,X2,X3,... vectors must match V.');
end
s = cell(size(y));
for i=1:length(s),
s{i} = 1 + (y{i}-x{i}(1))/(x{i}(xprodsiz(i))-x{i}(1))*(siz(i)-1);
end
end
for i=length(s):-1:1, ssiz{i} = size(s{i}); end
if ~isequal(ssiz{:}),
error('Y1,Y2,Y3, etc. must be the same size.');
end
% Check for out of range values of s and set to 1
sout = cell(size(s));
for i=1:length(s),
sout{i} = find((s{i}<.5) | (s{i}>=siz(i)+.5));
if ~isempty(sout{i}), s{i}(sout{i}) = ones(size(sout{i})); end
end
% Matrix element indexing
offset = cumprod([1 siz(1:end-1)]);
ndx = 1;
for i=1:length(s),
ndx = ndx + offset(i)*(round(s{i})-1);
end
% Now interpolate
F = v(ndx);
% Now set out of range values to NaN.
for i=1:length(sout),
if ~isempty(sout{i}), F(sout{i}) = NaN; end
end
%------------------------------------------------------
function [msg,x,v,y] = xnchk(x,v,y)
%XNCHK Check arguments to N-D data routines.
% [MSG,X,V,Y] = XNCHK(X,V,Y), checks the input aguments and
% returns either an error message in MSG or valid X,V, and Y
% data. X and Y are cell array lists that contain vectors or
% matrices.
if nargin~=3, error('Wrong number of input arguments.'); end
msg = [];
% Check to make sure the number of dimensions of V matches the
% number of x and y arguments.
if length(x) ==length(y) & ndims(v) ~= length(x),
msg = sprintf('V must be a %d-D array.',length(x));
return
end
siz = size(v);
isvec = zeros(size(x));
for i=length(x):-1:1,
xsiz{i} = size(x{i});
isvec(i) = isvector(x{i});
end
if ~isvector(v), % v is not a vector or scalar
% Convert x,y,z to row, column, and page matrices if necessary.
if all(isvec),
[x{1:end}] = ndgrid(x{:}); % Grid the vectors
for i=length(x):-1:1, xsiz{i} = size(x{i}); end
if ~isequal(siz,xsiz{:}),
msg = 'The lengths of X1,X2,X3,etc. must match the size of V.';
return
end
elseif any(isvec),
msg = 'X1,X2,X3, etc. must all be vectors or all be arrays.';
return
else
if ~isequal(siz,xsiz{:}),
msg = 'Arrays X1,X2,X3, etc. must be the same size as V.';
return
end
end
elseif isvector(v) % v is a vector
for i=length(x):-1:1, xlength(i) = length(x{i}); end
if any(~isvec),
msg = 'X1,X2,X3, etc. must be vectors when V is.';
return
elseif ~isequal(length(v),xlength{:}),
msg = 'X1,X2,X3, etc. must be the same length as V.';
return
end
end
isvec = zeros(size(y));
for i=length(y):-1:1,
ysiz{i} = size(y{i});
isvec(i) =isvector(y{i});
end
% If y1,y2,y3, etc. don't all have the same orientation, then
% build y1,y2,y3, etc. arrays.
if automesh(y{:})
[y{1:end}] = ndgrid(y{:});
elseif ~isequal(ysiz{:}),
msg = sprintf('Y1,Y2,Y3, etc. must all the same size or vectors of different orientations.');
end
function tf = isvector(x)
%ISVECTOR True if x has only one non-singleton dimension.
tf = length(x)==prod(size(x));
%----------------------------------------------------------
function F = splinen(x,v,xi)
%N-D spline interpolation
% Determine abscissa vectors
for i=1:length(x);
ind(1:length(x)) = {1}; ind{i} = ':';
x{i} = reshape(x{i}(ind{:}),1,size(x{i},i));
end
%
% Check for plaid data.
%
isplaid = 1;
for i=1:length(xi),
ind(1:length(xi)) = {1}; ind{i} = ':';
siz = size(xi{i}); siz(i) = 1;
xxi{i} = xi{i}(ind{:});
if ~isequal(repmat(xxi{i},siz),xi{i})
isplaid = 0;
end
xxi{i} = xxi{i}(:).'; % Make it a row
end
if ~isplaid
F = splncore(x,v,xi);
else
F = splncore(x,v,xxi,'gridded');
end
11. ФУНКЦИЯ unwrap (корректировка разрывов фазовых углов)
unwrap
Correct phase angles
Syntax
Q = unwrap(P)
Q = unwrap(P,tol)
Q = unwrap(P,[],dim)
Q = unwrap(P,tol,dim)
Description
Q = unwrap(P) corrects the radian phase angles in array P by adding multiples of ±2
when absolute jumps between consecutive array elements are greater than
radians. If P is a matrix, unwrap operates columnwise. If P is a multidimensional
array, unwrap operates on the first nonsingleton dimension.
Q = unwrap(P,tol) uses a jump tolerance tol instead of the default value, .
Q = unwrap(P,[],dim) unwraps along dim using the default tolerance.
Q = unwrap(P,tol,dim) uses a jump tolerance of tol.
Examples
Array P features smoothly increasing phase angles except for discontinuities at
elements (3,1) and (1,2).
P =
0 7.0686 1.5708 2.3562
0.1963 0.9817 1.7671 2.5525
6.6759 1.1781 1.9635 2.7489
0.5890 1.3744 2.1598 2.9452
The function Q = unwrap(P) eliminates these discontinuities.
Q =
0 0.7854 1.5708 2.3562
0.1963 0.9817 1.7671 2.5525
0.3927 1.1781 1.9635 2.7489
0.5890 1.3744 2.1598 2.9452
Limitations
The unwrap function detects branch cut crossings, but it can be fooled by sparse,
rapidly changing phase values.
See Also
abs, angle
Литература.
1. Анализ и обработка данных. Специальный справочник. Санкт-Петербург. Изд. Дом:
“Питер”, 2001 г.
2. Комплексирование геофизических методов.
3. Инструкция по электроразведке. Л. Недра.1984 г.
4. Справочник геофизика. Электроразведка. М. Недра.1991 г
5. Теория полей, применяемых в разведочной геофизике. Л.М. Альпин, Д.С.Даев, А.Д. Каринский. М.Недра. 1985 г.
6. Ваньян Л.Л. Электромагнитные зондирования. М. Научный мир.1997 г.
Не нашли, что искали? Воспользуйтесь поиском по сайту:
©2015 - 2024 stydopedia.ru Все материалы защищены законодательством РФ.
|