diff options
Diffstat (limited to 'macros/hilbert1.sci')
| -rw-r--r-- | macros/hilbert1.sci | 143 |
1 files changed, 105 insertions, 38 deletions
diff --git a/macros/hilbert1.sci b/macros/hilbert1.sci index fbcb136..476c00c 100644 --- a/macros/hilbert1.sci +++ b/macros/hilbert1.sci @@ -1,39 +1,106 @@ -function h= hilbert1(f, varargin) -//Analytic extension of real valued signal. -//Calling Sequence -// h= hilbert1(f) -// h= hilbert1(f,N) -// h= hilbert1(f,N,dim) -//Parameters -//f: real or complex valued scalar or vector -//N: The result will have length N -//dim : It analyses the input in this dimension -//Description -//h = hilbert1 (f) computes the extension of the real valued signal f to an analytic signal. If f is a matrix, the transformation is applied to each column. For N-D arrays, the transformation is applied to the first non-singleton dimension. -// -//real (h) contains the original signal f. imag (h) contains the Hilbert transform of f. -// -//hilbert1 (f, N) does the same using a length N Hilbert transform. The result will also have length N. -// -//hilbert1 (f, [], dim) or hilbert1 (f, N, dim) does the same along dimension dim. -//Examples -//## notice that the imaginary signal is phase-shifted 90 degrees -// t=linspace(0,10,256); -// z = hilbert1(sin(2*pi*0.5*t)); -// grid on; plot(t,real(z),';real;',t,imag(z),';imag;'); - -funcprot(0); -rhs= argn(2); -if(rhs<1 | rhs>3) - error("Wrong number of Input Arguments") -end - -select(rhs) - case 1 then - h= callOctave("hilbert", f); - case 2 then - h= callOctave("hilbert", f, varargin(1)); - case 3 then - h= callOctave("hilbert", f, varargin(1), varargin(2)); -end +/*Calling Sequence + h = hilbert1 (f) + h = hilbert1 (f, N) + h = hilbert1 (f, N, dim) +Description + Analytic extension of real valued signal. + h = hilbert (f) computes the extension of the real valued signal f to an analytic signal. + If f is a matrix, the transformation is applied to each column. + For N-D arrays, the transformation is applied to the first non-singleton dimension. + real (h) contains the original signal f. imag (h) contains the Hilbert transform of f. + hilbert (f, N) does the same using a length N Hilbert transform. The result will also have length N. + hilbert (f, [], dim) or hilbert (f, N, dim) does the same along dimension dim. +Dependencies + fft1, ifft1, ipermute +Example + //the magnitude of the hilbert transform eliminates the carrier + t=linspace(0,10,1024); + x=5*cos(0.2*t).*sin(100*t); + plot(t,x,t,abs(hilbert(x))); + */ +function f=hilbert1(f, N ,dim ) + // ------ PRE: initialization and dimension shifting --------- + nargin = argn(2); + if (nargin<1 || nargin>3) + error("Please enter valid number of inputs") + end + if ~isreal(f) + warning ('HILBERT: ignoring imaginary part of signal'); + f = real (f); + end + D=ndims(f); + select nargin + case 1 then + N=[]; + dim=[]; + case 2 then + dim=[] + end + // Dummy assignment. + order=1; + if isempty(dim) + dim=1; + if sum(size(f)>1)==1 + // We have a vector, find the dimension where it lives. + dim=find(size(f)>1); + end + else + if (length(dim)~=1 || ~or(type(dim)==[1 5 8])) + error('HILBERT: dim must be a scalar.'); + end + if modulo(dim,1)~=0 + error('HILBERT: dim must be an integer.'); + end + if (dim<1) || (dim>D) + error('HILBERT: dim must be in the range from 1 to %d.',D); + end + end + if (length(N)>1 || ~or(type(N)==[1 5 8])) + error('N must be a scalar.'); + elseif (~isempty(N) && modulo(N,1)~=0) + error('N must be an integer.'); + end + if dim>1 + order=[dim, 1:dim-1,dim+1:D]; + // Put the desired dimension first. + f=permute(f,order); + end + Ls=size(f,1); + // If N is empty it is set to be the length of the transform. + if isempty(N) + N=Ls; + end + // moduloember the exact size for later and modify it for the new length + permutedsize=size(f); + permutedsize(1)=N; + // Reshape f to a matrix. + f=resize_matrix(f,size(f,1),length(f)/size(f,1)); + W=size(f,2); + if ~isempty(N) + siz=size(f); + siz(1)=N; + f=resize_matrix(f,siz); + end + // ------- actual computation ----------------- + if N>2 + f=fft1(f); + if modulo(N,2)==0 + f=[f(1,:); + 2*f(2:N/2,:); + f(N/2+1,:); + zeros(N/2-1,W)]; + else + f=[f(1,:); + 2*f(2:(N+1)/2,:); + zeros((N-1)/2,W)]; + end + f=ifft1(f); + end + // ------- POST: Restoration of dimensions ------------ + // Restore the original, permuted shape. + f=matrix(f,permutedsize); + if dim>1 + // Undo the permutation. + f=ipermute(f,order); + end endfunction |