octave_packages/signal-1.1.3/hilbert.m

   1 ## Copyright (C) 2000 Paul Kienzle  <pkienzle@users.sf.net>
   2 ## Copyright (C) 2007 Peter L. Soendergaard
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## -*- texinfo -*-
  18 ## @deftypefn {Function File} {@var{h} =} hilbert (@var{f},@var{N},@var{dim})
  19 ## Analytic extension of real valued signal
  20 ##
  21 ## @code{@var{h}=hilbert(@var{f})} computes the extension of the real
  22 ## valued signal @var{f} to an analytic signal. If @var{f} is a matrix,
  23 ## the transformation is applied to each column. For N-D arrays,
  24 ## the transformation is applied to the first non-singleton dimension.
  25 ##
  26 ## @code{real(@var{h})} contains the original signal @var{f}.
  27 ## @code{imag(@var{h})} contains the Hilbert transform of @var{f}.
  28 ##
  29 ## @code{hilbert(@var{f},@var{N})} does the same using a length @var{N}
  30 ## Hilbert transform. The result will also have length @var{N}.
  31 ##
  32 ## @code{hilbert(@var{f},[],@var{dim})} or
  33 ## @code{hilbert(@var{f},@var{N},@var{dim})} does the same along
  34 ## dimension dim.
  35 ## @end deftypefn
  36
  37 function f=hilbert(f, N = [], dim = [])
  38
  39   % ------ PRE: initialization and dimension shifting ---------
  40
  41   if (nargin<1 || nargin>3)
  42     print_usage;
  43   end
  44
  45   if ~isreal(f)
  46     warning ('HILBERT: ignoring imaginary part of signal');
  47     f = real (f);
  48   end
  49
  50   D=ndims(f);
  51
  52   % Dummy assignment.
  53   order=1;
  54
  55   if isempty(dim)
  56     dim=1;
  57
  58     if sum(size(f)>1)==1
  59       % We have a vector, find the dimension where it lives.
  60       dim=find(size(f)>1);
  61     end
  62
  63   else
  64     if (numel(dim)~=1 || ~isnumeric(dim))
  65       error('HILBERT: dim must be a scalar.');
  66     end
  67     if rem(dim,1)~=0
  68       error('HILBERT: dim must be an integer.');
  69     end
  70     if (dim<1) || (dim>D)
  71       error('HILBERT: dim must be in the range from 1 to %d.',D);
  72     end
  73
  74   end
  75
  76   if (numel(N)>1 || ~isnumeric(N))
  77     error('N must be a scalar.');
  78   elseif (~isempty(N) && rem(N,1)~=0)
  79     error('N must be an integer.');
  80   end
  81
  82   if dim>1
  83     order=[dim, 1:dim-1,dim+1:D];
  84
  85     % Put the desired dimension first.
  86     f=permute(f,order);
  87
  88   end
  89
  90   Ls=size(f,1);
  91
  92   % If N is empty it is set to be the length of the transform.
  93   if isempty(N)
  94     N=Ls;
  95   end
  96
  97   % Remember the exact size for later and modify it for the new length
  98   permutedsize=size(f);
  99   permutedsize(1)=N;
 100
 101   % Reshape f to a matrix.
 102   f=reshape(f,size(f,1),numel(f)/size(f,1));
 103   W=size(f,2);
 104
 105   if ~isempty(N)
 106     f=postpad(f,N);
 107   end
 108
 109   % ------- actual computation -----------------
 110   if N>2
 111     f=fft(f);
 112
 113     if rem(N,2)==0
 114       f=[f(1,:);
 115          2*f(2:N/2,:);
 116          f(N/2+1,:);
 117          zeros(N/2-1,W)];
 118     else
 119       f=[f(1,:);
 120          2*f(2:(N+1)/2,:);
 121          zeros((N-1)/2,W)];
 122     end
 123
 124     f=ifft(f);
 125   end
 126
 127   % ------- POST: Restoration of dimensions ------------
 128
 129   % Restore the original, permuted shape.
 130   f=reshape(f,permutedsize);
 131
 132   if dim>1
 133     % Undo the permutation.
 134     f=ipermute(f,order);
 135   end
 136
 137 endfunction
 138
 139 %!demo
 140 %! % notice that the imaginary signal is phase-shifted 90 degrees
 141 %! t=linspace(0,10,256);
 142 %! z = hilbert(sin(2*pi*0.5*t));
 143 %! grid on; plot(t,real(z),';real;',t,imag(z),';imag;');
 144
 145 %!demo
 146 %! % the magnitude of the hilbert transform eliminates the carrier
 147 %! t=linspace(0,10,1024);
 148 %! x=5*cos(0.2*t).*sin(100*t);
 149 %! grid on; plot(t,x,'g;z;',t,abs(hilbert(x)),'b;|hilbert(z)|;');