function y = dst1( u )
% The discrete sine transform of type 1.
%
% This code also works for complex inputs. 
%
% To obtain the orthonormal eigenvectors of 
% 
% [ 2  -1             ]
% [ -1  2  -1         ]
% [                   ]
% [           -1  2 -1]
% [              -1  2]
% 
% try S = dst1( eye( N ) );  S = S*sqrt(2/(N+1)).
% 
% Alex Townsend, August 2014. Updated Sept 2014.

% Fold out to a FFT:
[n, m] = size(u);
% Be careful! DST has (n+1)(k+1) in its definition. Pad by a zero.  
Z = zeros(1, m);
% sin(kz) = (e^(ikz)-e^(-ikz))/2i 
% The term e^(-ikz)  =>  roots of unity in reverse order, flip!
y = fft( [ Z ; u ; Z ; -u(end:-1:1,:)]/2i );
% Take one of the halves: 
y = y( 2*n+2:-1:n+3, : );

end