1 files changed, 61 insertions, 29 deletions
diff --git a/macros/fht.sci b/macros/fht.sci
index f7e5c8e..9bd6b38 100644
--- a/macros/fht.sci
+++ b/macros/fht.sci
@@ -1,29 +1,61 @@
-function y=fht(d,n,dim)
-//The Function calculates the Fast Hartley Transform of real input.
-//Calling Sequence
-//M = fht (D)
-//M = fht (D, N)
-//M = fht (D, N, DIM)
-//Parameters 
-//Description
-//This function calculates the Fast Hartley transform of real input D. If D is a matrix, the Hartley transform is calculated along the columns by default.
-//Examples
-//fht(1:4)
-//ans =
-//   10   -4   -2   0  
-//This function is being called from Octave.
-funcprot(0);
-rhs=argn(2);
-if(rhs<1 | rhs>3)
-    error("Wrong number of input arguments.")
-end
-select(rhs)
-case 1 then
-    y=callOctave("fht",d)
-case 2 then
-    y=callOctave("fht",d,n)
-case 3 then
-    y=callOctave("fht",d,n,dim)
-end
-
-endfunction
+/* Description
+        This function calculates the Fast Hartley transform of real input D.
+        If D is a matrix, the Hartley transform is calculated along the columns by default.
+        If N is specified, the first N elements along the specified dimension are used for the transform.
+        If DIM is specified, the transform is calculated along the specified dimension.
+   Calling Sequence
+        M = fht (D)
+        M = fht (D, N)
+        M = fht (D, N, DIM)
+   Parameters 
+        D: Input data (real matrix or vector).
+        N: Number of elements of D to be used for the transform (optional).
+        DIM: Dimension along which the transform is to be computed (optional).
+   Examples
+        fht(1:4)
+        ans =
+       10   -4   -2   0  
+*/
+function M = fht(D, N, DIM)
+    funcprot(0);
+    rhs = argn(2);
+    if rhs < 1 | rhs > 3 then
+        error("Wrong number of input arguments.")
+    end
+    // The fht will be calculated along the first non-singleton dimension of the array i.e along the columns by default.
+    dimension = size(D);
+    nsdim = 1;
+    for i = 1:length(dimension)
+        if dimension(i) ~= 1 then
+            nsdim = i;
+            break;
+        end
+    end
+    // Process input arguments
+    select(rhs)
+    case 1 then
+        M = fft(D, -1, nsdim);
+    case 2 then
+        if isempty(N) then
+            n = size(D, nsdim);
+        else
+            n = N;
+        end
+        new_size = size(D);
+        new_size(nsdim) = n;
+        D = resize_matrix(D, new_size);
+        M = fft(D, -1, nsdim);
+    case 3 then
+        if isempty(N) then
+            n = size(D, DIM);
+        else
+            n = N;
+        end
+        new_size = size(D);
+        new_size(DIM) = n;
+        D = resize_matrix(D, new_size);
+        M = fft(D, -1, DIM);
+    end
+    // Return real part of the result
+    M = real(M) - imag(M);
+endfunction