001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 *  image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause 
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 
006 * Copyright (c) 2006-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 
007 * Visit http://imagingbook.com for additional details.
008 *******************************************************************************/
009
010package imagingbook.pub.dct;
011
012import imagingbook.lib.math.Matrix;
013
014/**
015 * This class calculates the 1D DFT using tabulated cosine values.
016 * This version is considerably faster than the original without tables.
017 * Other optimizations are possible.
018 * @author W. Burger
019 * @version 2015-07-29
020 */
021public class Dct1d {
022        
023        static final double CM0 = 1.0 / Math.sqrt(2);
024
025        final private double s;                 // common scale factor
026        final private double[] tmp;             // array to hold temporary data
027        final private int M;                    // size of the input vector
028        final private int N;                    // = 4M
029        
030        /**
031         * This table holds the cosine values cos(PI * (m(2u + 1)) / (2M)) = cos(j * PI / (2M))
032         * for all possible values of j = m (2u + 1), i.e., cosTable[j] =  cos[j * PI / (2M)].
033         * The table is of size N = 4M. To retrieve the correct cosine value for a given index 
034         * pair (m, u) use cosTable[(m * (2 * u + 1)) % (4 * M)].
035         */
036        final private double[] cosTable;
037        
038        public Dct1d(int M) {
039                this.M = M;
040                this.N = 4 * M;
041                this.s = Math.sqrt(2.0 / M); 
042                this.tmp = new double[M];
043                this.cosTable = makeCosineTable(M);
044        }
045        
046        private double[] makeCosineTable(int M) {
047                double[] table = new double[N];         // we need a table of size 4*M
048                for (int j = 0; j < N; j++) {           // j is equivalent to (m * (2 * u + 1)) % 4M
049                        double phi = j * Math.PI / (2 * M);
050                        table[j] = Math.cos(phi);
051                }
052                return table;
053        }
054
055        /**
056         * Performs the 1D forward discrete cosine transform.
057         * Destructively applies the forward Discrete Cosine Transform to 
058         * the argument vector.
059         * @param g the data to be transformed
060         */
061        public void DCT(double[] g) {
062                if (g.length != M)
063                        throw new IllegalArgumentException();
064                double[] G = tmp;
065                for (int m = 0; m < M; m++) {
066                        double cm = (m == 0) ? CM0 : 1.0;
067                        double sum = 0;
068                        for (int u = 0; u < M; u++) {
069                                sum += g[u] * cm * cosTable[(m * (2 * u + 1)) % N];
070                        }
071                        G[m] = s * sum;
072                }
073                System.arraycopy(G, 0, g, 0, M); // copy G -> g
074        }
075        
076        /**
077         * Performs the 1D inverse discrete cosine transform.
078         * Destructively applies the inverse Discrete Cosine Transform to 
079         * the argument vector.
080         * @param G the data to be transformed
081         */
082        public void iDCT(double[] G) {
083                if (G.length != M)
084                        throw new IllegalArgumentException();
085                double[] g = tmp;
086                for (int u = 0; u < M; u++) {
087                        double sum = 0;
088                        for (int m = 0; m < M; m++) {
089                                double cm = (m == 0) ? CM0 : 1.0;
090                                sum += cm * G[m] * cosTable[(m * (2 * u + 1)) % N];
091                        }
092                        g[u] = s * sum;
093                }
094                System.arraycopy(g, 0, G, 0, M); // copy g -> G
095        }
096        
097        
098        // test example
099        public static void main(String[] args) {
100                
101                double[] data = {1,2,3,4,5,3,0};
102                System.out.println("Original data:      " + Matrix.toString(data));
103                System.out.println();
104                
105                Dct1d dct = new Dct1d(data.length);
106                dct.DCT(data);
107                System.out.println("DCT of data:        " + Matrix.toString(data));
108                System.out.println();
109                
110                dct.iDCT(data);
111                System.out.println("Reconstructed data: " + Matrix.toString(data));
112                System.out.println();
113        }
114        
115//      Original data:      {1.000, 2.000, 3.000, 4.000, 5.000, 3.000, 0.000}
116//      DCT of data:        {6.803, -0.361, -3.728, 1.692, -0.888, -0.083, 0.167}
117//      Reconstructed data: {1.000, 2.000, 3.000, 4.000, 5.000, 3.000, -0.000}
118
119        
120}