Scilab / Octaveの比較 - アダマール変換

	  setdefaultlanguage('en_US')
	

	  // makematrix_hadamard.sci
	  // Copyright (C) 2008-2009 - INRIA - Michael Baudin
	  //
	  // This file must be used under the terms of the CeCILL.
	  // This source file is licensed as described in the file COPYING, which
	  // you should have received as part of this distribution.  The terms
	  // are also available at
	  // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

	  //
	  // makematrix_hadamard --
	  //   Returns the Hadamard matrix of size n.
	  //   n is expected to be a power of 2.
	  // Arguments, input
	  //  n : the size of the matrix to return
	  // Reference
	  //  http://en.wikipedia.org/wiki/Hadamard_matrix
	  //
	  function a = makematrix_hadamard ( n )
	    pow2n = log2(n);
	    if ( n - 2^pow2n <>0 ) then
	      errmsg = msprintf(gettext("%s: n is expected to be a power of 2."), "makematrix_hadamard");
	      error(errmsg)
	    end
	    if n==1 then
	      a = ones(1,1)
	    else
	      m = int(n/2)
	      b = makematrix_hadamard(m)
	      a = zeros(n,n)
	      a(1:m,1:m) = b(1:m,1:m)
	      a(1:m,m+1:n) = b(1:m,1:m)
	      a(m+1:n,1:m) = b(1:m,1:m)
	      a(m+1:n,m+1:n) = -b(1:m,1:m)
	    end
	  endfunction

	  // hadamard.sci
	  // Copyright (C) 2008-2009 - INRIA - Michael Baudin
	  //
	  // This file must be used under the terms of the CeCILL.
	  // This source file is licensed as described in the file COPYING, which
	  // you should have received as part of this distribution.  The terms
	  // are also available at
	  // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

	  //
	  // hadamard --
	  //   Returns the Hadamard matrix of size n.
	  //   n is expected to be a power of 2.
	  // Arguments, input
	  //  n : the size of the matrix to return
	  // Reference
	  //  http://en.wikipedia.org/wiki/Hadamard_matrix
	  //
	  function a = hadamard ( n )
	    a = makematrix_hadamard ( n );
	  endfunction
	

	  function hmatrix=hadamard(n)
	  h=[1,1;1,-1]
	  hmatrix=h;
	    for n = 1:1:(log2(n)+1)-2
	      // kron() - Kronecker product (.*.)
	      hmatrix = hmatrix.*.h;
	    end
	  endfunction
	

	  dec2bin(0:7)-'0'
	  (dec2bin(0:7)-'0')'
	  mod((dec2bin(0:7)-'0') * (dec2bin(0:7)-'0')', 2)
	  mod((dec2bin(0:7)-'0') * (dec2bin(0:7)-'0')', 2) * (-2) + 1
	

	  // Gray code
	  function list = graycode(N)
	    P=2^N;
	    list = dec2bin(bitxor(int((0:P-1)/2),0:P-1))';
	  endfunction
	  
	  graycode(3)
	  graycode(4)
	

	  % Gray code	  
	  function list = graycode(N)
	    P=2^N;
	    list = dec2bin(bitxor(bitshift(0:P-1,-1),0:P-1));
	  end
	  
	  graycode(3)
	  graycode(4)
	

	  // Length of Walsh (Hadamard) functions
	  N = 8;
	  hadamardMatrix=hadamard(N)

	  // Hadamard index
	  HadIdx = 0:N-1;
	  // Number of bits to represent the index
	  M = log2(N)+1;

	  //binHadIdx = fliplr(dec2bin(HadIdx,M))-'0';
	  //関数 flipdim(A,1) は flipud(A) と同じ
	  //関数 flipdim(A,2) は関数 fliplr(A) と同じ

	  // Bit reversing of the binary index
	  binHadIdx = strrev(dec2bin(HadIdx(:),M))
	  // Pre-allocate memory
	  binSeqVdx = zeros(N,M);
	  for k = 1:1:N
	      binSeqVdx(k,:) = strtod(strsplit(binHadIdx(k)))';
	  end

	  // Pre-allocate memory
	  binSeqIdx = zeros(N,M-1);
	  //for k = M:-1:2
	  //  binSeqIdx(:,k) = xor(binHadIdx(:,k),binHadIdx(:,k-1));
	  //end
	  for k = M:-1:2
	    binSeqIdx(:,k) = bitxor(binSeqVdx(:,k),binSeqVdx(:,k-1));
	  end

	  // Binary to integer sequency index
	  //SeqIdx = binSeqIdx*(2^((M-1:-1:0)'));
	  for k = 1:1:N
	    // Binary sequency index 
	    HadIdx(k) = bin2dec(strcat(string(binSeqIdx(k,:))));
	  end

	  // 1-based indexing
	  walshMatrix = hadamardMatrix(HadIdx+1,:)
	

	  pkg load signal

	  N = 8;  % Length of Walsh (Hadamard) functions
	  hadamardMatrix = hadamard(N)

	  % Hadamard index
	  HadIdx = 0:N-1;
	  % Number of bits to represent the index
	  M = log2(N)+1;

	  % Bit reversing of the binary index
	  binHadIdx = fliplr(dec2bin(HadIdx,M))-'0';
	  % Pre-allocate memory
	  binSeqIdx = zeros(N,M-1);
	  for k = M:-1:2
	    % Binary sequency index 
	    binSeqIdx(:,k) = xor(binHadIdx(:,k),binHadIdx(:,k-1));
	  end
	  % Binary to integer sequency index
	  SeqIdx = binSeqIdx*pow2((M-1:-1:0)');
	  % 1-based indexing
	  walshMatrix = hadamardMatrix(SeqIdx+1,:)
	

	  N = 8;
	  H = hadamard(N);
	  x = 8.*H(1,:) + 12.*H(3,:) + 18.*H(5,:) + 10.*H(8,:);
	  // fwht(x,8,'hadamard')
	  y = x * hadamardMatrix / N
	  // ifwht(y,8,'hadamard')
	  z = y * hadamardMatrix
	  norm(x-z)
	

	  % Fast Walsh-Hadamard transform
	  y1 = fwht(walshMatrix)
	  
	  N = 8;
	  H = hadamard(N);
	  x = 8.*H(1,:) + 12.*H(3,:) + 18.*H(5,:) + 10.*H(8,:);
	  y = fwht(x,N,'hadamard')
	  xHat = ifwht(y,N,'hadamard')
	  norm(x-xHat)
	

	  HadIdx
	

	  x = [19 -1 11 -9 -7 13 -15 5];
	  // y = fwht(x)
	  y = x * walshMatrix / 8
	

	  % Sequency (Walsh) Ordered Walsh-Hadamard Matrix
	  pkg load signal
	  N=8;
	  bitand(bitxor(bitshift(bitrevorder(0:N-1),1),bitrevorder(0:N-1)),(N-1))
	  bitand(bitrevorder(bitxor(bitshift(0:N-1,1),0:N-1)),(N-1))

	  % Dyadic (Paley) ordered Matrix
	  bitrevorder(0:N-1)

	  paleyMatrix = hadamardMatrix(bitrevorder(0:N-1)+1,:)
	  fwht(paleyMatrix,N,'dyadic')
	

	  x = [19 -1 11 -9 -7 13 -15 5];
	  y = fwht(x)
	  z = fwht(x,8,'dyadic')
	  p = x*paleyMatrix/8
	  porder = bitrevorder(0:N-1);
	  fwht(x(porder+1),N,'hadamard')
	  fwht(x,N,'hadamard')(porder+1)
	

	  function data=fhtnat(data)
	  N = length(data);
	  tmp = zeros(1,N);
	  level = 1;
	  for level=1:log2(N)
	    addbin = repmat([ones(1,2^(level-1)), zeros(1,2^(level-1))],1,2^(log2(N)-level));
	    addindex = find(addbin==1);
	    subindex = find(addbin==0);
	    tmp(addindex) = data(addindex) + data(subindex);
	    tmp(subindex) = data(addindex) - data(subindex);
	    data = tmp;
	  end
	  data = data/N;
	  endfunction
	  
	  function data=fhtseq(data)
	  N = length(data);
	  tmp = zeros(1,N);
	  addindex2 = find(repmat([1,0],1,N/2)==1);
	  subindex2 = find(repmat([1,0],1,N/2)==0);
	  for level=log2(N):-1:2
	    addbin = repmat([ones(1,2^(level-1)), zeros(1,2^(level-1))],1,2^(log2(N)-level));
	    addindex = find(addbin==1);
	    subindex = find(addbin==0);
	  
	    decbin = repmat([1,0,0,1],1,N/4);
	    addindexf = find(decbin==1);
	    subindexf = find(decbin==0);
	    
	    tmp(addindex) = data(addindex2) + data(subindex2);
	    tmp(subindex) = data(addindexf) - data(subindexf);
	    data = tmp;
	  end
	  tmp(addindex2) = data(addindex2) + data(subindex2);
	  tmp(subindex2) = data(addindex2) - data(subindex2);
	  data = tmp;
	  
	  data = data/N;
	  endfunction
	

	  pkg load signal;
	  x = [19 -1 11 -9 -7 13 -15 5];
	  y = fwht(x,length(x),'hadamard')
	  y = fhtnat(x)
	  y = fwht(x,length(x),'sequency')
	  y = fhtseq(x)
	

	  function data=fhtseqv2(data)
	  N = length(data);
	  tmp = zeros(1,N);

	  addindex2 = find(repmat([1,0],1,N/2)==1);
	  subindex2 = find(repmat([1,0],1,N/2)==0);
	  tmp(addindex2) = data(addindex2) + data(subindex2);
	  tmp(subindex2) = data(addindex2) - data(subindex2);
	  data = tmp;

	  addindexf = find(repmat([1,0,0,1],1,N/4)==1);
	  subindexf = find(repmat([1,0,0,1],1,N/4)==0);

	  for level=2:1:log2(N)
	    vecbin = repmat([ones(1,2^(level-1)), zeros(1,2^(level-1))],1,2^(log2(N)-level));
	    addindex = find(vecbin==1);
	    subindex = find(vecbin==0);
	    tmp(addindexf) = data(addindex) + data(subindex);
	    tmp(subindexf) = data(addindex) - data(subindex);
	    data = tmp;
	  end
	  data = data/N;
	  endfunction

	  pkg load signal;
	  x = [19 -1 11 -9 -7 13 -15 5];
	  y = fwht(x,length(x),'sequency')
	  y = fhtseqv2(x)
	

	  def FWHT(x):
	      """ Fast Walsh-Hadamard Transform
	      Based on mex function written by Chengbo Li@Rice Uni for his TVAL3 algorithm.
	      His code is according to the K.G. Beauchamp's book -- Applications of Walsh and Related Functions.
	      """
	      x = x.squeeze()
	      N = x.size
	      G = N/2 # Number of Groups
	      M = 2 # Number of Members in Each Group

	      # First stage
	      y = np.zeros((N/2,2))
	      y[:,0] = x[0::2] + x[1::2]
	      y[:,1] = x[0::2] - x[1::2]
	      x = y.copy()
	      # Second and further stage
	      for nStage in xrange(2,int(log(N,2))+1):
	          y = np.zeros((G/2,M*2))
	          y[0:G/2,0:M*2:4] = x[0:G:2,0:M:2] + x[1:G:2,0:M:2]
	          y[0:G/2,1:M*2:4] = x[0:G:2,0:M:2] - x[1:G:2,0:M:2]
	          y[0:G/2,2:M*2:4] = x[0:G:2,1:M:2] - x[1:G:2,1:M:2]
	          y[0:G/2,3:M*2:4] = x[0:G:2,1:M:2] + x[1:G:2,1:M:2]
	          x = y.copy()
	          G = G/2
	          M = M*2
	      x = y[0,:]
	      x = x.reshape((x.size,1))
	      return x/float(N)
	

	  /*=================================================================
	   *
	   * \file fWHtrans.cpp
	   *
	   *
	   * This code computes the (real) fast discrete Walsh-Hadamard transform with sequency order according to the K.G. Beauchamp's book -- Applications of Walsh and Related Functions.
	   *
	   *
	   * This file is written by Chengbo Li from Computational and Applied Mathematics Department of Rice University.
	   *
	   *
	   * This is a MEX-file for MATLAB.
	   *
	   * 02/15/2010
	   *
	   *=================================================================*/
	  #include <math.h>
	  #include "mex.h"
	  #include "matrix.h"
	  //#include <stdio.h>
	  //#include <stdlib.h>
	  // #include <malloc.h>
	  // #include <stack>


	  //! Matlab entry function
	  /*!
	   * \param nlhs number of left-hand-side output arguments
	   * \param plhs mxArray of output arguments
	   * \param nrhs number of right-hand-side input arguments
	   * \param prhs mxArray of input arguments
	   */
	  void mexFunction( int nlhs, mxArray *plhs[],
	          int nrhs, const mxArray*prhs[] )
	          
	  {
	      int p, nStage, L, clm;
	      int i, j, m, n, N, J, K, M;
	      double *v_in, *v_out, *v_ext, *v_final, *temp;
	      
	      /* Check for proper number of arguments */
	      if (nrhs != 1) {
	          mexErrMsgTxt("Only one input arguments required.");
	      }
	      else if (nlhs > 1) {
	          mexErrMsgTxt("Too many output arguments.");
	      }
	      
	      /* Get the size and pointers to input data. */
	      m  = mxGetM(prhs[0]);
	      n  = mxGetN(prhs[0]);
	      
	      
	      /* Make sure that both input and output vectors have the length with 2^p where p is some integer. */
	      p = ceil(log2(m));
	      N = pow(2, p);
	      
	      plhs[0] = mxCreateDoubleMatrix(N, n, mxREAL);
	      
	      v_final = mxGetPr(plhs[0]);
	      v_in = mxGetPr(prhs[0]);
	      
	      /* Extend the input vector if necessary. */
	      v_ext = (double*) mxCalloc(N, sizeof(double));
	      v_out = (double*) mxCalloc(N, sizeof(double));
	      
	      for (clm=0; clm<n; clm++) {
	          
	          /* C is row major while Matlab is column major */
	          for (j=0; j<m; j++){
	              v_ext[j] = (v_in+clm*m)[j];
	          }
	          for (j=m; j<N; j++){
	              v_ext[j] = 0;
	          }
	          
	          
	          for (i=0; i<N-1; i = i+2) {
	              v_ext[i] = v_ext[i] + v_ext[i+1];
	              v_ext[i+1] = v_ext[i] - v_ext[i+1] * 2;
	          }
	          L = 1;
	          
	          /* main loop */
	          for (nStage = 2; nStage<=p; ++nStage){
	              M = pow(2, L);
	              J = 0;
	              K = 0; // difference between Matlab and C++
	              while (K<N-1){
	                  for (j = J;j < J+M-1; j = j+2){
	                      /* sequency order  */
	                      v_out[K] = v_ext[j] + v_ext[j+M];
	                      v_out[K+1] = v_ext[j] - v_ext[j+M];
	                      v_out[K+2] = v_ext[j+1] - v_ext[j+1+M];
	                      v_out[K+3] = v_ext[j+1] + v_ext[j+1+M];
	                      
	                      K = K+4;
	                  }
	                  
	                  J = J+2*M;
	              }
	              
	              temp = v_ext;
	              v_ext = v_out;
	              v_out = temp;
	              
	              L = L+1;
	          }
	          
	          /* Perform scaling of coefficients. */
	          for ( i =0; i<N; ++i){
	              (v_final+clm*N)[i] = v_ext[i]/N;
	          }
	      }
	      
	      /* Set free the memory. */
	      mxFree(v_out);
	      mxFree(v_ext);
	      
	      return;
	  }
	

	  # -*- coding: utf-8 -*-
	  """
	  Created on Thu Jun 11 15:21:43 2015
	  Fast Walsh-Hadamard Transform with Sequency Order
	  Author: Ding Luo@Fraunhofer IOSB
	  """
	  from math import log
	  import numpy as np
	  import GrayCode
	  from time import clock

	  def get_sequency_list(inputArray):
	      """ Sort input 1D array into sequency order
	      Utilizes gray code generation from a Python recipe from Internet.
	      """
	      length = inputArray.size
	      bitlength = int(log(length,2))
	      # Gray Code
	      graycodes=GrayCode.GrayCode(bitlength)
	      # Bitreverse of gray code
	      bitreverse = [int(graycodes[i][::-1],2) for i in range(length)]
	      
	      outputArray = inputArray.copy()
	      outputArray[bitreverse] = inputArray[:]

	      return outputArray

	  def SFWHT(X):
	      """ 'Slower' Fast Walsh-Hadamard Transform
	      Step#1 Get sequency-ordered input
	      Step#2 Perform Hadamard Transform
	      """
	      x=get_sequency_list(X)
	      N = x.size
	      M = int(log(N,2))
	      out = x.copy()
	      for m in range(M):
	          outtemp = out.copy()         
	          step = 2**m
	          numCalc = 2**m
	          for g in range(0,N,2*step): # number of groups
	              for c in range(numCalc):
	                  index = g + c
	                  out[index] = outtemp[index] + outtemp[index+step]
	                  out[index+step] = outtemp[index] - outtemp[index+step]
	      return out/float(N)


	  def FWHT(x):
	      """ Fast Walsh-Hadamard Transform
	      Based on mex function written by Chengbo Li@Rice Uni for his TVAL3 algorithm.
	      His code is according to the K.G. Beauchamp's book -- Applications of Walsh and Related Functions.
	      """
	      x = x.squeeze()
	      N = x.size
	      G = N // 2 # Number of Groups
	      M = 2 # Number of Members in Each Group

	      # First stage
	      y = np.zeros((G,2))
	      y[:,0] = x[0::2] + x[1::2]
	      y[:,1] = x[0::2] - x[1::2]
	      x = y.copy()
	      # Second and further stage
	      for nStage in range(2,int(log(N,2))+1):
	          y = np.zeros((G//2,M*2))
	          y[0:G//2,0:M*2:4] = x[0:G:2,0:M:2] + x[1:G:2,0:M:2]
	          y[0:G//2,1:M*2:4] = x[0:G:2,0:M:2] - x[1:G:2,0:M:2]
	          y[0:G//2,2:M*2:4] = x[0:G:2,1:M:2] - x[1:G:2,1:M:2]
	          y[0:G//2,3:M*2:4] = x[0:G:2,1:M:2] + x[1:G:2,1:M:2]
	          x = y.copy()
	          G = int(G/2)
	          M = M*2
	      x = y[0,:]
	      x = x.reshape((x.size,1))
	      return x/float(N)

	     
	  if __name__ == "__main__":
	      x = np.random.random(1024**2)
	      t1 = clock()
	      y1 = SFWHT(x)
	      t1 = clock() - t1
	      print(t1)
	      t2 = clock()
	      y2 = FWHT(x)
	      t2 = clock() - t2
	      print(t2)
	

	  def isOdd(integer):
	      #assert isinstance(integer, int)
	      return integer % 2 == 1

	  def isEven(integer):
	      #assert isinstance(integer, int)
	      return integer % 2 == 0

	  def _list_to_string(li):
	      return ''.join(map(str, li))

	  class GrayCode(object):
	      def __init__(self, nbits):
	          self._nbits = nbits
	          self._grayCode = []
	          self.__generate()

	      def __getitem__(self, i):
	          return self._grayCode[i]

	      def __str__(self):
	          return str(self._grayCode)

	      __repr__ = __str__

	      def __iter__(self):            
	          return self._grayCode.__iter__()

	      def __generate(self):
	          li = [0 for i in range(self._nbits)]
	          self._grayCode.append(_list_to_string(li))

	          for term in range(2, (1<<self._nbits)+1):
	              if isOdd(term):                
	                  for i in range(-1,-(self._nbits),-1):
	                      if li[i]==1:                        
	                          li[i-1]=li[i-1]^1                        
	                          break
	                      
	              if isEven(term):
	                  li[-1]=li[-1]^1

	              self._grayCode.append(_list_to_string(li))

	  class GrayCodeIterator(object):
	      def __init__(self, nbits):
	          self._nbits = nbits

	      def __iter__(self):
	          li = [0 for i in range(self._nbits)]
	          yield _list_to_string(li)

	          for term in range(2, (1<<self._nbits)+1):
	              if isOdd(term):                
	                  for i in range(-1,-(self._nbits),-1):
	                      if li[i]==1:                        
	                          li[i-1]=li[i-1]^1                        
	                          break
	                      
	              if isEven(term):
	                  li[-1]=li[-1]^1

	              yield _list_to_string(li)


	  if __name__=='__main__':
	      d = 4
	      codes=GrayCode(d)
	      print ('%d digits gray codes:', d)
	      print (codes)
	      print ('Using Iterator:')
	      #for c in GrayCode(20):  
	      #    print c
	      for c in GrayCodeIterator(20):
	          print (c)
	

	  import FWHT as fwht
	  import GrayCode as gc
	  import numpy as np
	  import importlib
	  importlib.reload(fwht)
	  importlib.reload(gc)
	  
	  a = np.array([19, -1, 11, -9, -7, 13, -15, 5])
	  
	  # Secuency ordered
	  print("Secuency ordered")
	  print(fwht.SFWHT(a))
	  print(fwht.FWHT(a).T)
	  
	  from scipy.linalg import hadamard
	  print(hadamard(4, dtype=complex))
	  print(hadamard(8))
	  # Hadamard ordered
	  print(np.dot(a,hadamard(8))/8)
	  
	  # // division error example
	  print("matrix indice")
	  print(a[0:8//2])
	  print(a[0:8/2])

	  print(np.zeros((8//2,2)))
	  np.zeros(8//2,2)
	

	  Secuency ordered
	  [  2.   3.   0.   4.   0.   0.  10.   0.]
	  [[  2.   3.   0.   4.   0.   0.  10.   0.]]
	  [[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]
	   [ 1.+0.j -1.-0.j  1.+0.j -1.-0.j]
	   [ 1.+0.j  1.+0.j -1.-0.j -1.-0.j]
	   [ 1.+0.j -1.-0.j -1.-0.j  1.+0.j]]
	  [[ 1  1  1  1  1  1  1  1]
	   [ 1 -1  1 -1  1 -1  1 -1]
	   [ 1  1 -1 -1  1  1 -1 -1]
	   [ 1 -1 -1  1  1 -1 -1  1]
	   [ 1  1  1  1 -1 -1 -1 -1]
	   [ 1 -1  1 -1 -1  1 -1  1]
	   [ 1  1 -1 -1 -1 -1  1  1]
	   [ 1 -1 -1  1 -1  1  1 -1]]
	  [  2.   0.   4.   0.   3.  10.   0.   0.]
	  matrix indice
	  [19 -1 11 -9]

	  ---------------------------------------------------------------------------
	  TypeError Traceback (most recent call last)
	  <ipython-input-194-830b68925cc9> in <module>()
	       22 print("matrix indice")
	       23 print(a[0:8//2])
	  ---> 24 print(a[0:8/2])

	  TypeError: slice indices must be integers or None or have an __index__ method
	

	  print(np.zeros((8//2,2)))
	  np.zeros(8//2,2)
	  
	  np.zeros(8/2,2)
	

	  [[ 0.  0.]
	   [ 0.  0.]
	   [ 0.  0.]
	   [ 0.  0.]]
	  
	  ---------------------------------------------------------------------------
	  TypeError Traceback (most recent call last)
	  <ipython-input-205-d3ad594fc4dc> in <module>()
	        1 print(np.zeros((8//2,2)))
	  ----> 2 np.zeros(8//2,2)
	  
	  TypeError: data type not understood
	  
	  ---------------------------------------------------------------------------
	  TypeError Traceback (most recent call last)
	  <ipython-input-206-805c1441d66f> in <module>()
	  ----> 1 np.zeros(8/2,2)
	  
	  TypeError: 'float' object cannot be interpreted as an integer
	

	  atomsInstall("IPCV")

	  N=256;
	  hadamardMatrix=hadamard(N);

	  function wmat=had2wal(hmat,msize)
	  HadIdx = 0:msize-1;
	  M = log2(msize)+1;
	  binHadIdx = strrev(dec2bin(HadIdx(:),M));
	  binSeqVdx = zeros(msize,M);
	  for k = 1:1:msize
	      binSeqVdx(k,:) = strtod(strsplit(binHadIdx(k)))';
	  end
	  binSeqIdx = zeros(msize,M-1);
	  for k = M:-1:2
	      binSeqIdx(:,k) = bitxor(binSeqVdx(:,k),binSeqVdx(:,k-1));
	  end
	  for k = 1:1:msize
	      HadIdx(k) = bin2dec(strcat(string(binSeqIdx(k,:))));
	  end
	  wmat = hmat(HadIdx+1,:);
	  endfunction
	  
	  walshMatrix = had2wal(hadamardMatrix,N);
	  imwrite(hadamardMatrix,"hmat.bmp");
	  imwrite(walshMatrix,"wmat.bmp");
	  
	  img_orig = double(imread("BRIDGE.bmp"));
	  imshow(uint8(img_orig));
	  surf(imdct(img_orig))
	  surf(imdct(img_orig)(1:50,1:50))
	

	  img_orig = double(imread("Text.bmp"));
	  imshow(uint8(img_orig));
	  img_dct = dct2(img_orig);
	  surf(img_dct);
	  surf(img_dct(1:50,1:50));
	

	  function imat = fwht2d(imat,hmat)
	  N = sqrt(length(imat))
	  tmat = zeros(N,N);
	  for i = 1:N
	    tmat(i,:) = imat(i,:) * hmat / N;
	  end
	  imat = zeros(N,N);
	  for j = 1:N
	    imat(:,j) = hmat * tmat(:,j)
	  end
	  endfunction
	  
	  a2d = matrix(1:64,8,8)'

	  // Sequency order
	  fwht2d(a2d,walshMatrix)
	  walshMatrix*a2d*walshMatrix / 8
	

	  // 通常は意味のない行為
	  img_dct = imdct(img_orig);
	  img_wal = fwht2d(img_dct,walshMatrix);
	  imshow(uint8( (img_wal-min(img_wal))/(max(img_wal)-min(img_wal))*255 ));
	  
	  img_mat = fwht2d(img_orig,walshMatrix);
	  img_ict = imidct(img_mat);
	  imshow(uint8( (img_ict-min(img_ict))/(max(img_ict)-min(img_ict))*255 ));