1 files changed, 65 insertions, 30 deletions
diff --git a/macros/spectral_xdf.sci b/macros/spectral_xdf.sci
index f0d457b..6081c71 100644
--- a/macros/spectral_xdf.sci
+++ b/macros/spectral_xdf.sci
@@ -1,31 +1,66 @@
-function y= spectral_xdf(x, varargin)
-// Return the spectral density estimator given a data vector X, window name WIN, and bandwidth, B.
-//Calling Sequence
-//spectral_xdf(X)
-//spectral_xdf(X, WIN)
-//spectral_xdf(X, WIN, B)
-//Parameters 
-//X: Data Vector
-//WIN: Window names
-//B: Bandwidth 
-//Description
-//Return the spectral density estimator given a data vector X, window name WIN, and bandwidth, B.
-//The window name, e.g., "triangle" or "rectangle" is used to search for a function called 'WIN_lw'.
-//If WIN is omitted, the triangle window is used.
-//If B is omitted, '1 / sqrt (length (X))' is used.
-funcprot(0);
-rhs= argn(2);
-if(rhs<1 | rhs>3)
-error("Wrong number of input arguments")
-end
-
-select(rhs)
-	case 1 then
-		y= callOctave("spectral_xdf", x);
-	case 2 then
-		y= callOctave("spectral_xdf", x , varargin(1));
-	case 3 then
-		y= callOctave("spectral_xdf", x , varargin(1), varargin(2));
-
-end
+/* Description
+        Return the spectral density estimator given a data vector x, window name win, and bandwidth, b.
+        The window name, e.g., "triangle" or "rectangle" is used to search for a function called win_sw.
+        If win is omitted, the triangle window is used.
+        If b is omitted, 1 / sqrt (length (x)) is used.
+    Calling Sequence
+        spectral_xdf (x)
+        spectral_xdf (x, win)
+        spectral_xdf (x, win, b)
+    Parameters
+        x : Data vector
+        win : the window name . Default "triangle" is used .
+        b : Bandwidth . Default value 1/sqrt(length(x))
+Dependencies:fft1 ifft1 */
+function sde = spectral_xdf (x, win, b)
+    // check x is a vector or not
+    if ~isvector(x) 
+        error("spectral_xdf :  x must a data vector")
+    end
+    // available windows - rectangle
+    function c = rectangle_sw (n, b)
+          c = zeros (n, 1);
+          c(1) = 2 / b + 1;
+          l = (2:n)' - 1;
+          l = 2 * %pi * l / n;
+          c(2:n) = sin ((2/b + 1) * l / 2) ./ sin (l / 2);
+    endfunction
+    // triangle
+    function retval = triangle_sw (n, b)
+          retval = zeros (n,1);
+          retval(1) = 1 / b;
+          l = (2:n)' - 1;
+          l = 2 * %pi * l / n;
+          retval(2:n) = b * (sin (l / (2*b)) ./ sin (l / 2)).^2;
+    endfunction
+    // main function
+    nargin = argn(2)
+    if (nargin < 1)
+      error("invalid inputs");
+    end
+    xr = max(size (x) );
+    if (size (x, 2) > 1)
+      x = x';
+    end
+    if (nargin < 3)
+      b = 1 / ceil (sqrt (xr));
+    end
+    if (nargin == 1)
+      w = triangle_sw (xr, b);
+    elseif (~ (type(win)== 10))
+      error ("spectral_xdf: WIN must be a string");
+    elseif (~strcmp (win , "triangle" ) ) 
+        w = triangle_sw (xr , b);
+    elseif (~strcmp (win , "rectangle" ) )     
+        w = rectangle_sw ( xr , b);
+    else
+      error("Invalid window or this window is not supported yet");
+    end
+    x = x - sum (x) / xr;
+    sde = (abs (fft1 (x)) / xr).^2;
+    sde = real (ifft1 (fft1 (sde) .* fft1 (w)));
+    zer = zeros(xr,1);
+    sde = [zer sde];
+    sde(:, 1) = (0 : xr-1)' / xr;
+    sde = clean(sde);
 endfunction