Tanggal :September 22, 2020

# py29. Discrete Cosine Transform in Python

### For comparison, we use the dct and idct functions from scipy.fftpack module with aliases of sdct, and sidct.

`# ex29.pyfrom __future__ import print_function, divisionimport numpy as npfrom scipy.fftpack import dct as sdct, idct as sidctfrom numpy.fft import fft, ifftimport matplotlib.pyplot as pltdef dct(X):    x = np.concatenate([X,np.flipud(X)])    N = 2*len(X)    f0 = 1/np.sqrt(2*N)    f1 = 1/np.sqrt(N)    Y = fft(x)    kin = np.arange(N)    Mul = np.cos(np.pi/N*kin)-1j*np.sin(np.pi/N*kin)    Y *= Mul    Y *= f1    Y *= f0/f1    return Y[:N/2].realdef idct(X):    N = len(X)    f0 = 1/np.sqrt(4*N)    f1 = 1/np.sqrt(2*N)    kin = np.arange(N)    FTR = X*np.cos(.5*np.pi/N*kin)    FTI = X*np.sin(.5*np.pi/N*kin)    FTR /= f1    FTI /= f1    FTR *= f1/f0    FTI *= f1/f0    FTR1 = np.concatenate((FTR,np.zeros(1),np.flipud(FTR[1:])))    FTI1 = np.concatenate((FTI,np.zeros(1),-1*np.flipud(FTI[1:])))    y = ifft(FTR1+1j*FTI1)    return y[:N].realif __name__ == '__main__':    x = np.array(range(4)+range(2,5)+range(10,1,-1))    y = dct(x)    y1 = sdct(x.astype('float'), norm = 'ortho')    z = idct(y)    z1 = sidct(y1, norm = 'ortho')    err1 = (z-x).std()    err2 = (z1-x).std()    print ('err1 =',err1)    print ('err2 =',err2)    plt.plot(x,'ro--', markerfacecolor='w')    plt.plot(y,'bo--', markerfacecolor='w')    plt.xlabel('n')    plt.ylabel('x,y')    plt.title('x (red), y = dct(x) (blue)')    plt.xlim(-.5,15.5)    plt.show()# err1 = 1.43874580946e-15# err2 = 4.52626554495e-16`

Share