1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
|
"""
Real spectrum tranforms (DCT, DST, MDCT)
"""
__all__ = ['dct', 'idct']
import numpy as np
from scipy.fftpack import _fftpack
from scipy.fftpack.basic import _datacopied
import atexit
atexit.register(_fftpack.destroy_ddct1_cache)
atexit.register(_fftpack.destroy_ddct2_cache)
atexit.register(_fftpack.destroy_dct1_cache)
atexit.register(_fftpack.destroy_dct2_cache)
def dct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=0):
"""
Return the Discrete Cosine Transform of arbitrary type sequence x.
Parameters
----------
x : array_like
The input array.
type : {1, 2, 3}, optional
Type of the DCT (see Notes). Default type is 2.
n : int, optional
Length of the transform.
axis : int, optional
Axis over which to compute the transform.
norm : {None, 'ortho'}, optional
Normalization mode (see Notes). Default is None.
overwrite_x : bool, optional
If True the contents of x can be destroyed. (default=False)
Returns
-------
y : ndarray of real
The transformed input array.
See Also
--------
idct
Notes
-----
For a single dimension array ``x``, ``dct(x, norm='ortho')`` is equal to
MATLAB ``dct(x)``.
There are theoretically 8 types of the DCT, only the first 3 types are
implemented in scipy. 'The' DCT generally refers to DCT type 2, and 'the'
Inverse DCT generally refers to DCT type 3.
type I
~~~~~~
There are several definitions of the DCT-I; we use the following
(for ``norm=None``)::
N-2
y[k] = x[0] + (-1)**k x[N-1] + 2 * sum x[n]*cos(pi*k*n/(N-1))
n=1
Only None is supported as normalization mode for DCT-I. Note also that the
DCT-I is only supported for input size > 1
type II
~~~~~~~
There are several definitions of the DCT-II; we use the following
(for ``norm=None``)::
N-1
y[k] = 2* sum x[n]*cos(pi*k*(2n+1)/(2*N)), 0 <= k < N.
n=0
If ``norm='ortho'``, ``y[k]`` is multiplied by a scaling factor `f`::
f = sqrt(1/(4*N)) if k = 0,
f = sqrt(1/(2*N)) otherwise.
Which makes the corresponding matrix of coefficients orthonormal
(``OO' = Id``).
type III
~~~~~~~~
There are several definitions, we use the following
(for ``norm=None``)::
N-1
y[k] = x[0] + 2 * sum x[n]*cos(pi*(k+0.5)*n/N), 0 <= k < N.
n=1
or, for ``norm='ortho'`` and 0 <= k < N::
N-1
y[k] = x[0] / sqrt(N) + sqrt(1/N) * sum x[n]*cos(pi*(k+0.5)*n/N)
n=1
The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up
to a factor `2N`. The orthonormalized DCT-III is exactly the inverse of
the orthonormalized DCT-II.
References
----------
http://en.wikipedia.org/wiki/Discrete_cosine_transform
'A Fast Cosine Transform in One and Two Dimensions', by J. Makhoul, `IEEE
Transactions on acoustics, speech and signal processing` vol. 28(1),
pp. 27-34, http://dx.doi.org/10.1109/TASSP.1980.1163351 (1980).
"""
if type == 1 and norm is not None:
raise NotImplementedError(
"Orthonormalization not yet supported for DCT-I")
return _dct(x, type, n, axis, normalize=norm, overwrite_x=overwrite_x)
def idct(x, type=2, n=None, axis=-1, norm=None, overwrite_x=0):
"""
Return the Inverse Discrete Cosine Transform of an arbitrary type sequence.
Parameters
----------
x : array_like
The input array.
type : {1, 2, 3}, optional
Type of the DCT (see Notes). Default type is 2.
n : int, optional
Length of the transform.
axis : int, optional
Axis over which to compute the transform.
norm : {None, 'ortho'}, optional
Normalization mode (see Notes). Default is None.
overwrite_x : bool, optional
If True the contents of x can be destroyed. (default=False)
Returns
-------
y : ndarray of real
The transformed input array.
See Also
--------
dct
Notes
-----
For a single dimension array `x`, ``idct(x, norm='ortho')`` is equal to
MATLAB ``idct(x)``.
'The' IDCT is the IDCT of type 2, which is the same as DCT of type 3.
IDCT of type 1 is the DCT of type 1, IDCT of type 2 is the DCT of type
3, and IDCT of type 3 is the DCT of type 2. For the definition of these
types, see `dct`.
"""
if type == 1 and norm is not None:
raise NotImplementedError(
"Orthonormalization not yet supported for IDCT-I")
# Inverse/forward type table
_TP = {1:1, 2:3, 3:2}
return _dct(x, _TP[type], n, axis, normalize=norm, overwrite_x=overwrite_x)
def _dct(x, type, n=None, axis=-1, overwrite_x=0, normalize=None):
"""
Return Discrete Cosine Transform of arbitrary type sequence x.
Parameters
----------
x : array-like
input array.
n : int, optional
Length of the transform.
axis : int, optional
Axis along which the dct is computed. (default=-1)
overwrite_x : bool, optional
If True the contents of x can be destroyed. (default=False)
Returns
-------
z : real ndarray
"""
tmp = np.asarray(x)
if not np.isrealobj(tmp):
raise TypeError("1st argument must be real sequence")
if n is None:
n = tmp.shape[axis]
else:
raise NotImplemented("Padding/truncating not yet implemented")
if tmp.dtype == np.double:
if type == 1:
f = _fftpack.ddct1
elif type == 2:
f = _fftpack.ddct2
elif type == 3:
f = _fftpack.ddct3
else:
raise ValueError("Type %d not understood" % type)
elif tmp.dtype == np.float32:
if type == 1:
f = _fftpack.dct1
elif type == 2:
f = _fftpack.dct2
elif type == 3:
f = _fftpack.dct3
else:
raise ValueError("Type %d not understood" % type)
else:
raise ValueError("dtype %s not supported" % tmp.dtype)
if normalize:
if normalize == "ortho":
nm = 1
else:
raise ValueError("Unknown normalize mode %s" % normalize)
else:
nm = 0
if type == 1 and n < 2:
raise ValueError("DCT-I is not defined for size < 2")
overwrite_x = overwrite_x or _datacopied(tmp, x)
if axis == -1 or axis == len(tmp.shape) - 1:
return f(tmp, n, nm, overwrite_x)
#else:
# raise NotImplementedError("Axis arg not yet implemented")
tmp = np.swapaxes(tmp, axis, -1)
tmp = f(tmp, n, nm, overwrite_x)
return np.swapaxes(tmp, axis, -1)
|
