function y = my_fftw_dst(x)
%
% Y = DST(X) returns the discrete Sine transform of X. The ouput
%     matches that FFTW's DST-I. If X is a matrix, each column 
%     is transformed.
%

if min(size(x)) == 1
   x = x(:);
end
[n m] = size(x);
yy = zeros(2*(n+1),m); 
yy(2:n+1,:) = x; 
yy(n+3:2*(n+1),:) = -flipud(x); % odd extension
yy = fft(yy); 
% recover sine terms 
y = -imag(yy(2:n+1,:));



 



