function Y = lp_dst1(X)
% LP_DST1  Discrete sine transform of type I along the columns of X.
%
% Y(k,:) = sqrt(2/(n+1)) * sum_j X(j,:) sin(pi*j*k/(n+1)).  The scaling makes
% the transform orthogonal, so lp_dst1(lp_dst1(X)) is X, and it is computed
% with one complex FFT of the odd extension of X, of length 2(n+1).
n = size(X, 1);
m = size(X, 2);
Z = zeros(2*(n+1), m);
Z(2:n+1, :) = X;
Z(n+3:2*(n+1), :) = -flipud(X);
W = fft(Z);
Y = -imag(W(2:n+1, :)) * sqrt(1/(2*(n+1)));
end