Xcorr.java

/*
 * $Id: Xcorr.java,v 1.19 2008/07/17 07:30:03 koga Exp $
 * 
 * Copyright (C) 2004 Koga Laboratory. All rights reserved.
 */
package org.mklab.tool.signal;

import java.util.List;

import org.mklab.nfc.matrix.DoubleMatrix;
import org.mklab.tool.matrix.Makecolv;


/**
 * 相関係数を求めるクラスです。
 * 
 * <p> Cross-correlation function
 * 
 * @author koga
 * @version $Revision: 1.19 $
 * @see org.mklab.tool.matrix.Cov
 * @see org.mklab.tool.signal.Xcov
 */
public class Xcorr {

  /**
   * もし<code>x</code>がベクトルなら、自己相関係数を返します。
   * 
   * <p> もし<code>x</code>が<code>m</code>×<code>n</code>の行列なら、 <code>n^2</code>列が<code>x</code>の列の全ての組み合わせについての 相互相関係数を含む行列を返します。
   * 
   * @param x_ データ
   * @return 自己相関係数 cross-correlation
   */
  public static DoubleMatrix xcorr(DoubleMatrix x_) {
    DoubleMatrix x = x_;

    DoubleMatrix X;
    if (Math.min(x.getColumnSize(), x.getRowSize()) == 1) {
      X = Makecolv.makecolv(x).appendRight(Makecolv.makecolv(x));
    } else {
      X = x;
    }

    int xm = X.getRowSize();
    int xn = X.getColumnSize();
    int m = 2 * xm - 1;
    DoubleMatrix c = x.createZero(m, xn * xn);

    for (int i = 1; i <= xn; i++) {
      DoubleMatrix tmp = X.getSubMatrix(xm, 1, i, i).conjugate();
      tmp.setElement(m, 0);
      for (int j = 1; j <= i; j++) {
        List<DoubleMatrix> temp = Filter.filter(X.getSubMatrix(1, X.getRowSize(), j, j), x.createUnit(1), tmp);
        DoubleMatrix cc = temp.get(0);
        int k1 = (i - 1) * xn + j;
        int k2 = (j - 1) * xn + i;
        c.setSubMatrix(1, c.getRowSize(), k2, k2, cc.transpose());
        if (i != j) {
          c.setSubMatrix(1, c.getRowSize(), k1, k1, c.getSubMatrix(m, 1, k2, k2).conjugate());
        }
      }
    }

    c = c.getSubMatrix(1, c.getRowSize(), 1, 1);

    if (x_.getRowSize() == 1) {
      c = c.transpose();
    }

    return c;

  }

  /**
   * <code>x</code>と<code>y</code>の<code>(2*m-1)</code>点の相互相関係数を返します。 ただし、<code>x</code>と<code>y</code>は<code>m</code>点のベクトルです。
   * 
   * @param x_ データ1
   * @param y_ データ2
   * @return 相互相関係数 (cross-correlation)
   */
  public static DoubleMatrix xcorr(DoubleMatrix x_, DoubleMatrix y_) {
    DoubleMatrix x = x_;
    DoubleMatrix y = y_;

    if (Math.min(x.getColumnSize(), x.getRowSize()) != 1 && Math.min(y.getColumnSize(), y.getRowSize()) != 1) {
      throw new IllegalArgumentException(Messages.getString("Xcorr.0")); //$NON-NLS-1$
    }

    int nx = x.length();
    int ny = y.length();
    if (nx > ny) {
      y.setElement(nx, 0);
    } else if (nx < ny) {
      x.setElement(ny, 0);
    }
    DoubleMatrix X = Makecolv.makecolv(x).appendRight(Makecolv.makecolv(y));

    int xm = X.getRowSize();
    int xn = X.getColumnSize();
    int m = 2 * xm - 1;
    DoubleMatrix c = x.createZero(m, xn * xn);

    for (int i = 1; i <= xn; i++) {
      DoubleMatrix tmp = X.getSubMatrix(xm, 1, i, i).conjugate();
      tmp.setElement(m, 0);
      for (int j = 1; j <= i; j++) {
        List<DoubleMatrix> temp = Filter.filter(X.getSubMatrix(1, X.getRowSize(), j, j), x.createUnit(1), tmp);
        DoubleMatrix cc = temp.get(0);
        int k1 = (i - 1) * xn + j;
        int k2 = (j - 1) * xn + i;
        c.setSubMatrix(1, c.getRowSize(), k2, k2, cc.transpose());
        if (i != j) {
          c.setSubMatrix(1, c.getRowSize(), k1, k1, c.getSubMatrix(m, 1, k2, k2).conjugate());
        }
      }
    }

    c = c.getSubMatrix(1, c.getRowSize(), 1, 1);

    if (x_.getRowSize() == 1) {
      c = c.transpose();
    }

    return c;
  }

}