|
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 disrete cosine transform of x. If n is given, then |
|
|
19 ## x is padded or trimmed to length n before computing the transform. |
|
|
20 ## If x is a matrix, compute the transform along the columns of the |
|
|
21 ## the matrix. The transform is faster if x is real-valued and even |
|
|
22 ## length. |
|
|
23 ## |
|
|
24 ## The discrete cosine transform X of x can be defined as follows: |
|
|
25 ## |
|
|
26 ## N-1 |
|
|
27 ## X[k] = w(k) sum x[n] cos (pi (2n-1) k / 2N ), k = 0, ..., N-1 |
|
|
28 ## n=0 |
|
|
29 ## |
|
|
30 ## with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1. There |
|
|
31 ## are other definitions with different scaling of X[k], but this form |
|
|
32 ## is common in image processing. |
|
|
33 ## |
|
|
34 ## See also: idct, dct2, idct2, dctmtx |
|
|
35 |
|
|
36 ## From Discrete Cosine Transform notes by Brian Evans at UT Austin, |
|
|
37 ## http://www.ece.utexas.edu/~bevans/courses/ee381k/lectures/09_DCT/lecture9/ |
|
|
38 ## the discrete cosine transform of x at k is as follows: |
|
|
39 ## |
|
|
40 ## N-1 |
|
|
41 ## X[k] = sum 2 x[n] cos (pi (2n-1) k / 2N ) |
|
|
42 ## n=0 |
|
|
43 ## |
|
|
44 ## which can be computed using: |
|
|
45 ## |
|
|
46 ## y = [ x ; flipud (x) ] |
|
|
47 ## Y = fft(y) |
|
|
48 ## X = exp( -j pi [0:N-1] / 2N ) .* Y |
|
|
49 ## |
|
|
50 ## or for real, even length x |
|
|
51 ## |
|
|
52 ## y = [ even(x) ; flipud(odd(x)) ] |
|
|
53 ## Y = fft(y) |
|
|
54 ## X = 2 real { exp( -j pi [0:N-1] / 2N ) .* Y } |
|
|
55 ## |
|
|
56 ## Scaling the result by w(k)/2 will give us the desired output. |
|
|
57 |
|
|
58 ## Author: Paul Kienzle |
|
|
59 ## 2001-02-08 |
|
|
60 ## * initial release |
|
|
61 function y = dct (x, n) |
|
|
62 |
|
|
63 if (nargin < 1 || nargin > 2) |
|
|
64 usage ("y = dct(x [, n])"); |
|
|
65 endif |
|
|
66 |
|
|
67 realx = isreal(x); |
|
|
68 transpose = (rows (x) == 1); |
|
|
69 |
|
|
70 if transpose, x = x (:); endif |
|
|
71 [nr, nc] = size (x); |
|
|
72 if nargin == 1 |
|
|
73 n = nr; |
|
|
74 elseif n > nr |
|
|
75 x = [ x ; zeros(n-nr,nc) ]; |
|
|
76 elseif n < nr |
|
|
77 x (nr-n+1 : n, :) = []; |
|
|
78 endif |
|
|
79 |
|
|
80 if n == 1 |
|
|
81 w = 1/2; |
|
|
82 else |
|
|
83 w = [ sqrt(1/4/n); sqrt(1/2/n)*exp((-1i*pi/2/n)*[1:n-1]') ] * ones (1, nc); |
|
|
84 endif |
|
|
85 if ( realx && rem (n, 2) == 0 ) |
|
|
86 y = fft([ x(1:2:n,:) ; x(n:-2:1,:) ]); |
|
|
87 y = 2 * real( w .* y ); |
|
|
88 else |
|
|
89 y = fft ([ x ; flipud (x) ]); |
|
|
90 y = w .* y (1:n, :); |
|
|
91 if (realx) y = real (y); endif |
|
|
92 endif |
|
|
93 if transpose, y = y.'; endif |
|
|
94 |
|
|
95 endfunction |
