#!/usr/bin/perl
######################
# Testing the DCT matrices from Ian's blur.cgi program.
# Jim on April 7
######################
use strict;
use warnings;
use Math::Trig;
use Math::Matrix;

# --- from Ian's blur.cgi ----------------------------------------------

# Generate a DCT matrix for dimensions (rows,columns). 
# This is defined in the Wikipedia as DCT-II.
sub get_dctm {
  my $n = shift;
  my $pi_div_n = pi/$n;
  my $dctm = new Math::Matrix ([]);
  for my $row (0..$n-1){
    for my $col (0..$n-1){
      $dctm->[$row][$col] = cos ($pi_div_n * $row * ($col + 0.5));
    }
  }
  return $dctm;
}

# Generate the inverse of the DCT. 
# This is defined in the Wikipedia as DCT-III multiplied by 2/n.
sub get_idctm {
  my $n = shift;
  my $pi_div_n = pi/$n;
  my $one_div_n = 1/$n;
  my $idctm = new Math::Matrix ([]);

  # Jim: Changed column index to start at 1.
  #      Replaced pi/$n with $pi_div_n
  for my $row (0..$n-1){
    for my $col (1..$n-1){
      $idctm->[$row][$col] = (2/$n) * cos ($pi_div_n * $col * ($row + 0.5));
    }
  }

  # Jim: The k=0 term is different; left column of matrix is just 0.5 * 2/$n
  for my $j (0..$n-1){
    $idctm->[$j][0] = $one_div_n;
  }
  
  return $idctm;
}

# -- tests ---------------------------------------------------------------

#use Test::More qw(no_plan);

my $n = 4;
my $dct  = get_dctm($n,$n);
my $idct = get_idctm($n,$n);

print " DCT with n=$n \n";
$dct->print;
print "\n";
print " inverse DCT \n";
$idct->print;
print "\n";
print " product of the two \n";
$dct->multiply($idct)->print;


# see http://en.wikipedia.org/wiki/Discrete_cosine_transform
