Merge pull request #17 from avinashlalotra/masterHEADmaster

Abinash's Work

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