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) 
示例测试结果如下：
 
 
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
 
运行结果如下：