[Scipy-tickets] [SciPy] #733: Implement dct

SciPy scipy-tickets@scipy....
Thu Jan 15 10:02:07 CST 2009


#733: Implement dct
---------------------------+------------------------------------------------
 Reporter:  cdavid         |        Owner:  cdavid  
     Type:  defect         |       Status:  assigned
 Priority:  high           |    Milestone:  0.8.0   
Component:  scipy.fftpack  |      Version:          
 Severity:  normal         |   Resolution:          
 Keywords:                 |  
---------------------------+------------------------------------------------
Comment (by pbrod):

 One possible implementation of dct and idct:
 {{{
 import numpy as np
 def dct(x,n=None):
     """
     Discrete Cosine Transform

                       N-1
            y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
                       n=0

     Examples
     --------
     >>> import numpy as np
     >>> x = np.arange(5)
     >>> np.abs(x-idct(dct(x)))<1e-14
     array([ True,  True,  True,  True,  True], dtype=bool)
     >>> np.abs(x-dct(idct(x)))<1e-14
     array([ True,  True,  True,  True,  True], dtype=bool)

     Reference
     ---------
     http://en.wikipedia.org/wiki/Discrete_cosine_transform
     http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
     """
     fft = np.fft.fft
     x = np.atleast_1d(x)

     if n is None:
         n = x.shape[-1]

     if x.shape[-1]<n:
         n_shape = x.shape[:-1] + (n-x.shape[-1],)
         xx = np.hstack((x,np.zeros(n_shape)))
     else:
         xx = x[...,:n]

     real_x = np.all(np.isreal(xx))
     if (real_x and (np.remainder(n,2) == 0)):
         xp = 2 * fft(np.hstack( (xx[...,::2], xx[...,::-2]) ))
     else:
         xp = fft(np.hstack((xx, xx[...,::-1])))
         xp = xp[...,:n]

     w = np.exp(-1j * np.arange(n) * np.pi/(2*n))

     y = xp*w

     if real_x:
         return y.real
     else:
         return y

 def idct(x,n=None):
     """
     Inverse Discrete Cosine Transform

                        N-1
            x[k] = 1/N sum w[n]*y[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
                        n=0

            w(0) = 1/2
            w(n) = 1 for n>0

     Examples
     --------
     >>> import numpy as np
     >>> x = np.arange(5)
     >>> np.abs(x-idct(dct(x)))<1e-14
     array([ True,  True,  True,  True,  True], dtype=bool)
     >>> np.abs(x-dct(idct(x)))<1e-14
     array([ True,  True,  True,  True,  True], dtype=bool)

     Reference
     ---------
     http://en.wikipedia.org/wiki/Discrete_cosine_transform
     http://users.ece.utexas.edu/~bevans/courses/ee381k/lectures/
     """

     ifft = np.fft.ifft
     x = np.atleast_1d(x)

     if n is None:
         n = x.shape[-1]

     w = np.exp(1j * np.arange(n) * np.pi/(2*n))

     if x.shape[-1]<n:
         n_shape = x.shape[:-1] + (n-x.shape[-1],)
         xx = np.hstack((x,np.zeros(n_shape)))*w
     else:
         xx = x[...,:n]*w

     real_x = np.all(np.isreal(x))
     if (real_x and (np.remainder(n,2) == 0)):
         xx[...,0] = xx[...,0]*0.5
         yp = ifft(xx)
         y  = np.zeros(xx.shape,dtype=complex)
         y[...,::2] = yp[...,:n/2]
         y[...,::-2] = yp[...,n/2::]
     else:
         yp = ifft(np.hstack((xx, np.zeros_like(xx[...,0]),
 np.conj(xx[...,:0:-1]))))
         y = yp[...,:n]

     if real_x:
         return y.real
     else:
         return y
 }}}

-- 
Ticket URL: <http://scipy.org/scipy/scipy/ticket/733#comment:4>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list