Add a useful package (from Source forge) for octave

[CreaPhase.git] / octave_packages / signal-1.1.3 / hilbert.m

diff --git a/octave_packages/signal-1.1.3/hilbert.m b/octave_packages/signal-1.1.3/hilbert.m
new file mode 100644 (file)
index 0000000..92fe16e
--- /dev/null
+++ b/octave_packages/signal-1.1.3/hilbert.m
@@ -0,0 +1,149 @@
+## Copyright (C) 2000 Paul Kienzle  <pkienzle@users.sf.net>
+## Copyright (C) 2007 Peter L. Soendergaard 
+##
+## This program is free software; you can redistribute it and/or modify it under
+## the terms of the GNU General Public License as published by the Free Software
+## Foundation; either version 3 of the License, or (at your option) any later
+## version.
+##
+## This program is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+## details.
+##
+## You should have received a copy of the GNU General Public License along with
+## this program; if not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{h} =} hilbert (@var{f},@var{N},@var{dim})
+## Analytic extension of real valued signal
+##
+## @code{@var{h}=hilbert(@var{f})} computes the extension of the real
+## valued signal @var{f} to an analytic signal. If @var{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.
+##
+## @code{real(@var{h})} contains the original signal @var{f}.
+## @code{imag(@var{h})} contains the Hilbert transform of @var{f}.
+##
+## @code{hilbert(@var{f},@var{N})} does the same using a length @var{N}
+## Hilbert transform. The result will also have length @var{N}.
+##
+## @code{hilbert(@var{f},[],@var{dim})} or
+## @code{hilbert(@var{f},@var{N},@var{dim})} does the same along
+## dimension dim.
+## @end deftypefn
+
+function f=hilbert(f, N = [], dim = [])
+  
+  % ------ PRE: initialization and dimension shifting ---------
+   
+  if (nargin<1 || nargin>3)
+    print_usage;
+  end
+
+  if ~isreal(f)
+    warning ('HILBERT: ignoring imaginary part of signal');
+    f = real (f);
+  end
+  
+  D=ndims(f);
+  
+  % 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 (numel(dim)~=1 || ~isnumeric(dim))
+      error('HILBERT: dim must be a scalar.');
+    end
+    if rem(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 (numel(N)>1 || ~isnumeric(N))
+    error('N must be a scalar.');
+  elseif (~isempty(N) && rem(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
+  
+  % Remember the exact size for later and modify it for the new length
+  permutedsize=size(f);
+  permutedsize(1)=N;
+  
+  % Reshape f to a matrix.
+  f=reshape(f,size(f,1),numel(f)/size(f,1));
+  W=size(f,2);
+  
+  if ~isempty(N)
+    f=postpad(f,N);
+  end
+  
+  % ------- actual computation ----------------- 
+  if N>2
+    f=fft(f);
+    
+    if rem(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=ifft(f);
+  end
+  
+  % ------- POST: Restoration of dimensions ------------
+  
+  % Restore the original, permuted shape.
+  f=reshape(f,permutedsize);
+  
+  if dim>1
+    % Undo the permutation.
+    f=ipermute(f,order);
+  end
+  
+endfunction
+
+%!demo
+%! % notice that the imaginary signal is phase-shifted 90 degrees
+%! t=linspace(0,10,256);
+%! z = hilbert(sin(2*pi*0.5*t));
+%! grid on; plot(t,real(z),';real;',t,imag(z),';imag;');
+
+%!demo
+%! % the magnitude of the hilbert transform eliminates the carrier
+%! t=linspace(0,10,1024);
+%! x=5*cos(0.2*t).*sin(100*t);
+%! grid on; plot(t,x,'g;z;',t,abs(hilbert(x)),'b;|hilbert(z)|;');