function b = dst(a)
%%% compute the discrete sine transform.
%%% no error checking is computed.  a had better be a column vector or a
%%% matrix, ideally a power of 2 in length!!!


[m,n]=size(a);

y=zeros(2*(n+1),m);
y(2:n+1,:)=a;
y(n+3:2*(n+1),:)=-flipud(a);
yy=fft(y);
b=yy(2:n+1,:)/(-2*sqrt(-1));

if isreal(a), b = real(b); end