Сделай Сам Свою Работу на 5

Необходимые программы-функции





 

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 Все материалы защищены законодательством РФ.