python 实现maltab中离散正弦变换dst及其反变换idst

DFT、DCT、DST数学公式及原理部分可参考该博客。
笔者在将某matlab程序代码重新用python实现过程中，发现scipy.fft.dst的运行结果和matlab中的dst函数运行结果不一致（如下图所示）。前者主要有type 和axis两个参数，改变参数也无法得到和matlab中相同结果，便用python重写了matlab中的dst函数及其反变换idst函数。欢迎大家批评指正！

• matlab 中dst 源码

function b=dst(a,n)
%DST    Discrete sine transform.
%
%       Y = DST(X) returns the discrete sine transform of X.
%       The vector Y is the same size as X and contains the
%       discrete sine transform coefficients.
%
%       Y = DST(X,N) pads or truncates the vector X to length N
%       before transforming.
%
%       If X is a matrix, the DST operation is applied to each
%       column. This transform can be inverted using IDST.
%
%       See also: FFT, IFFT, IDST.

%       A. Nordmark 10-25-94.
%       Copyright 1994-2016 The MathWorks, Inc.


narginchk(1,2);

if min(size(a))==1
    if size(a,2)>1
        do_trans = 1;
    else
        do_trans = 0;
    end
    a = a(:);
else
    do_trans = 0;
end
if nargin==1
  n = size(a,1);
end
m = size(a,2);

% Pad or truncate a if necessary
if size(a,1)<n
  aa = zeros(n,m);
  aa(1:size(a,1),:) = a;
else
  aa = a(1:n,:);
end

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

if isreal(a), b = real(b); end
if do_trans, b = b.'; end

% LocalWords:  Nordmark
%%%%%%%%%%%%%%%%%%%%%%%%%%%
function b=idst(a,n)
%IDST   Inverse discrete sine transform.
%
%       X = IDST(Y) inverts the DST transform, returning the
%       original vector if Y was obtained using Y = DST(X).
%
%       X = IDST(Y,N) pads or truncates the vector Y to length N
%       before transforming.
%
%       If Y is a matrix, the IDST operation is applied to
%       each column.
%
%       See also: FFT, IFFT, DST.

%       A. Nordmark 10-25-94.
%       Copyright 1994-2001 The MathWorks, Inc.

if nargin==1
  if min(size(a))==1
    n=length(a);
  else
    n=size(a,1);
  end
end

nn=n+1;
b=2/nn*dst(a,n);

示例测试结果如下：

• dst的python实现

def DST(a,*args):
    shape = a.shape
    if len(shape) == 1:
        if shape[0] > 1:
            do_trans = 1
        else:
            do_trans = 0
        a = np.ravel(a)
        a = a.reshape(-1,1)
    else:
        do_trans = 0

    shape = a.shape

    if args:
        n = args[0]
    else:
        n = shape[0]

    m = a.shape[1]

    if shape[0] < n:
        aa = np.zeros((n,m))
        aa[0:shape[0],:] = a
    else:
        aa = a[0:n,:]

    y = np.zeros((2 * (n + 1),m))
    y[1:n+1,:] = aa

    y[n+2:2*(n+1),:] = -np.flipud(aa)
    yy = np.fft.fft(y,axis=0)

    b = yy[1:n+1,:] / (-2j)

    if np.isreal(a.any()):
        b = np.real(b)
    if do_trans:
        b = b.T

    return b

def IDST(a,*args):
    if args:
        n = args[0]
    else:
        n = a.shape[0]

    nn = n + 1
    b = 2/nn * DST(a,n)
    return b

运行结果如下：