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.
016 * Explicit (slow) version for that does not use pre-calculated cosine tables.
017 * 
018 * @author W. Burger
019 * @version 2014-04-13 (changed method signatures to operate destructively on the input arrays)
020 *
021 */
022public class Dct1d_Slow {
023        
024        static final double CM0 = 1.0 / Math.sqrt(2);
025
026        final private double[] tmp;
027        final private int M;
028        
029        public Dct1d_Slow(int M) {
030                this.M = M;
031                this.tmp = new double[M];
032        }
033        
034        /**
035         * Destructively applies the forward Discrete Cosine Transform to argument vector.
036         * @param g the data to be transformed
037         */
038        public void DCT(double[] g) {
039                if (g.length != M)
040                        throw new IllegalArgumentException();
041                final double s = Math.sqrt(2.0 / M); 
042                double[] G = tmp;
043                for (int m = 0; m < M; m++) {
044                        double cm = (m == 0) ? CM0 : 1.0;
045                        double sum = 0;
046                        for (int u = 0; u < M; u++) {
047                                double Phi = Math.PI * m * (2 * u + 1) / (2 * M);
048                                sum += g[u] * cm * Math.cos(Phi);
049                        }
050                        G[m] = s * sum;
051                }
052                System.arraycopy(G, 0, g, 0, M); // copy G -> g
053        }
054        
055        /**
056         * Destructively applies the inverse Discrete Cosine Transform to argument vector.
057         * @param G the data to be transformed
058         */
059        public void iDCT(double[] G) {
060                if (G.length != M)
061                        throw new IllegalArgumentException();
062                final double s = Math.sqrt(2.0 / M); //common scale factor
063                double[] g = tmp;
064                for (int u = 0; u < M; u++) {
065                        double sum = 0;
066                        for (int m = 0; m < M; m++) {
067                                double cm = (m == 0) ? CM0 : 1.0;
068                                double Phi = Math.PI * m * (2 * u + 1) / (2 * M);
069                                sum += G[m] * cm * Math.cos(Phi);
070                        }
071                        g[u] = s * sum;
072                }
073                System.arraycopy(g, 0, G, 0, M); // copy g -> G
074        }
075        
076        
077        // test example
078        public static void main(String[] args) {
079                
080                double[] data = {1,2,3,4,5,3,0};
081                System.out.println("Original data:      " + Matrix.toString(data));
082                
083                Dct1d_Slow dct = new Dct1d_Slow(data.length);
084                dct.DCT(data);
085                System.out.println("DCT of data:        " + Matrix.toString(data));
086                
087                dct.iDCT(data);
088                System.out.println("Reconstructed data: " + Matrix.toString(data));
089        }
090        
091//      Original data:      {1.000, 2.000, 3.000, 4.000, 5.000, 3.000, 0.000}
092//      DCT of data:        {6.803, -0.361, -3.728, 1.692, -0.888, -0.083, 0.167}
093//      Reconstructed data: {1.000, 2.000, 3.000, 4.000, 5.000, 3.000, -0.000}
094
095}