function b=dst1(a,n)
%DST1   Unitary Discrete Sine Transform Type-I.
%       (Adapted from MATLAB's internal dst.m function)
%
%       Y = DST1(X) returns the discrete sine transform of X.
%       The vector Y is the same size as X and contains the
%       discrete sine transform coefficients.
%
%       Y = DST1(X,N) pads or truncates the vector X to length N
%       before transforming.
%
%       If X is a matrix, the DST operation is applied to each
%       column. This transform can be inverted by applying the DST again.
%
%       See also: FFT, IFFT, IDST.
%       A. Nordmark 10-25-94. 
%       Copyright 1994-2003 The MathWorks, Inc.
%       $Revision: 1.8.4.1 $  $Date: 2003/11/18 03:11:20 $


error(nargchk(1,2,nargin));

if min(size(a))==1
    if size(a,2)>1
        do_trans = 1;
    else
        do_trans = 0;
    end
    a = a(:);
else
    do_trans = 0;
end
if nargin==1,
  n = size(a,1);
end
m = size(a,2);

% Pad or truncate a if necessary
if size(a,1)<n,
  aa = zeros(n,m);
  aa(1:size(a,1),:) = a;
else
  aa = a(1:n,:);
end

% Check for powers of 2
if mod(log2(n+1),1)
    warning(sprintf('DST1: n=%d. Pick a power of 2 for most efficiency\n',n+1));
end
if mod(log2(m+1),1)
    warning(sprintf('DST1: m=%d. Pick a power of 2 for most efficiency\n',m+1));
end

% OLD CODE BELOW
% y=zeros(2*(n+1),m+1); % m+1 very important, in order to get back power of 2.
% y(2:n+1,1:end-1)=aa;
% y(n+3:2*(n+1),1:end-1)=-flipud(aa);
y = [[zeros(1,m); aa; zeros(1,m); -flipud(aa)],zeros(2*n+2,1)];
yy=fft(y);
b=imag(yy(2:n+1,1:end-1));

if isreal(a), b = real(b); end
if do_trans, b = b.'; end

% Set this transform to be unitary!
b = -sqrt(1/2/(n+1))*b;
