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