|
0
|
1 ## Copyright (C) 2001 Paul Kienzle |
|
|
2 ## |
|
|
3 ## This program is free software; you can redistribute it and/or modify |
|
|
4 ## it under the terms of the GNU General Public License as published by |
|
|
5 ## the Free Software Foundation; either version 2 of the License, or |
|
|
6 ## (at your option) any later version. |
|
|
7 ## |
|
|
8 ## This program is distributed in the hope that it will be useful, |
|
|
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|
|
11 ## GNU General Public License for more details. |
|
|
12 ## |
|
|
13 ## You should have received a copy of the GNU General Public License |
|
|
14 ## along with this program; if not, write to the Free Software |
|
|
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA |
|
|
16 |
|
|
17 ## y = dct (x, n) |
|
|
18 ## Computes the inverse discrete cosine transform of x. If n is |
|
|
19 ## given, then x is padded or trimmed to length n before computing |
|
|
20 ## the transform. If x is a matrix, compute the transform along the |
|
|
21 ## columns of the the matrix. The transform is faster if x is |
|
|
22 ## real-valued and even length. |
|
|
23 ## |
|
|
24 ## The inverse discrete cosine transform x of X can be defined as follows: |
|
|
25 ## |
|
|
26 ## N-1 |
|
|
27 ## x[n] = sum w(k) X[k] cos (pi (2n-1) k / 2N ), k = 0, ..., N-1 |
|
|
28 ## k=0 |
|
|
29 ## |
|
|
30 ## with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1 |
|
|
31 ## |
|
|
32 ## See also: idct, dct2, idct2, dctmtx |
|
|
33 |
|
|
34 ## Author: Paul Kienzle |
|
|
35 ## 2001-02-08 |
|
|
36 ## * initial release |
|
|
37 function y = idct (x, n) |
|
|
38 |
|
|
39 if (nargin < 1 || nargin > 2) |
|
|
40 usage ("y = dct(x [, n])"); |
|
|
41 endif |
|
|
42 |
|
|
43 realx = isreal(x); |
|
|
44 transpose = (rows (x) == 1); |
|
|
45 |
|
|
46 if transpose, x = x (:); endif |
|
|
47 [nr, nc] = size (x); |
|
|
48 if nargin == 1 |
|
|
49 n = nr; |
|
|
50 elseif n > nr |
|
|
51 x = [ x ; zeros(n-nr,nc) ]; |
|
|
52 elseif n < nr |
|
|
53 x (n-nr+1 : n, :) = []; |
|
|
54 endif |
|
|
55 |
|
|
56 if ( realx && rem (n, 2) == 0 ) |
|
|
57 w = [ sqrt(n/4); sqrt(n/2)*exp((1i*pi/2/n)*[1:n-1]') ] * ones (1, nc); |
|
|
58 y = ifft (w .* x); |
|
|
59 y([1:2:n, n:-2:1], :) = 2*real(y); |
|
|
60 elseif n == 1 |
|
|
61 y = x; |
|
|
62 else |
|
|
63 ## reverse the steps of dct using inverse operations |
|
|
64 ## 1. undo post-fft scaling |
|
|
65 w = [ sqrt(4*n); sqrt(2*n)*exp((1i*pi/2/n)*[1:n-1]') ] * ones (1, nc); |
|
|
66 y = x.*w; |
|
|
67 |
|
|
68 ## 2. reconstruct fft result and invert it |
|
|
69 w = exp(-1i*pi*[n-1:-1:1]'/n) * ones(1,nc); |
|
|
70 y = ifft ( [ y ; zeros(1,nc); y(n:-1:2,:).*w ] ); |
|
|
71 |
|
|
72 ## 3. keep only the original data; toss the reversed copy |
|
|
73 y = y(1:n, :); |
|
|
74 if (realx) y = real (y); endif |
|
|
75 endif |
|
|
76 if transpose, y = y.'; endif |
|
|
77 |
|
|
78 endfunction |
