1 ## Author: Paul Kienzle <pkienzle@users.sf.net> (2006)
2 ## This program is granted to the public domain.
5 ## @deftypefn {Function File} {@var{y} =} dst (@var{x})
6 ## @deftypefnx {Function File} {@var{y} =} dst (@var{x}, @var{n})
7 ## Computes the type I discrete sine transform of @var{x}. If @var{n} is given,
8 ## then @var{x} is padded or trimmed to length @var{n} before computing the transform.
9 ## If @var{x} is a matrix, compute the transform along the columns of the
12 ## The discrete sine transform X of x can be defined as follows:
16 ## X[k] = sum x[n] sin (pi n k / (N+1) ), k = 1, ..., N
23 function y = dst (x, n)
25 if (nargin < 1 || nargin > 2)
29 transpose = (rows (x) == 1);
30 if transpose, x = x (:); endif
36 x = [ x ; zeros(n-nr,nc) ];
38 x (nr-n+1 : n, :) = [];
41 y = fft ([ zeros(1,nc); x ; zeros(1,nc); -flipud(x) ])/-2j;
43 if isreal(x), y = real (y); endif
45 ## Compare directly against the slow transform
47 # w = pi*[1:n]'/(n+1);
48 # for k = 1:n, y2(k) = sum(x(:).*sin(k*w)); end
51 if transpose, y = y.'; endif
56 %! x = log(linspace(0.1,1,32));
58 %! assert(y(3), sum(x.*sin(3*pi*[1:32]/33)), 100*eps)