Created
March 4, 2017 14:12
-
-
Save wyhasany/5cacb533a976a6eb748554532772e1f2 to your computer and use it in GitHub Desktop.
Simple cross correlation implementation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/****************************************************************************** | |
* Compilation: javac Complex.java | |
* Execution: java Complex | |
* | |
* Data type for complex numbers. | |
* | |
* The data type is "immutable" so once you create and initialize | |
* a Complex object, you cannot change it. The "final" keyword | |
* when declaring re and im enforces this rule, making it a | |
* compile-time error to change the .re or .im fields after | |
* they've been initialized. | |
* | |
* % java Complex | |
* a = 5.0 + 6.0i | |
* b = -3.0 + 4.0i | |
* Re(a) = 5.0 | |
* Im(a) = 6.0 | |
* b + a = 2.0 + 10.0i | |
* a - b = 8.0 + 2.0i | |
* a * b = -39.0 + 2.0i | |
* b * a = -39.0 + 2.0i | |
* a / b = 0.36 - 1.52i | |
* (a / b) * b = 5.0 + 6.0i | |
* conj(a) = 5.0 - 6.0i | |
* |a| = 7.810249675906654 | |
* tan(a) = -6.685231390246571E-6 + 1.0000103108981198i | |
* | |
******************************************************************************/ | |
/** | |
* The {@code Complex} class represents a complex number. | |
* Complex numbers are immutable: their values cannot be changed after they | |
* are created. | |
* It includes methods for addition, subtraction, multiplication, division, | |
* conjugation, and other common functions on complex numbers. | |
* <p> | |
* For additional documentation, see <a href="http://algs4.cs.princeton.edu/99scientific">Section 9.9</a> of | |
* <i>Algorithms, 4th Edition</i> by Robert Sedgewick and Kevin Wayne. | |
* | |
* @author Robert Sedgewick | |
* @author Kevin Wayne | |
*/ | |
public class Complex { | |
private final double re; // the real part | |
private final double im; // the imaginary part | |
/** | |
* Initializes a complex number from the specified real and imaginary parts. | |
* | |
* @param real the real part | |
* @param imag the imaginary part | |
*/ | |
public Complex(double real, double imag) { | |
re = real; | |
im = imag; | |
} | |
/** | |
* Returns a string representation of this complex number. | |
* | |
* @return a string representation of this complex number, | |
* of the form 34 - 56i. | |
*/ | |
public String toString() { | |
if (im == 0) return re + ""; | |
if (re == 0) return im + "i"; | |
if (im < 0) return re + " - " + (-im) + "i"; | |
return re + " + " + im + "i"; | |
} | |
/** | |
* Returns the absolute value of this complex number. | |
* This quantity is also known as the <em>modulus</em> or <em>magnitude</em>. | |
* | |
* @return the absolute value of this complex number | |
*/ | |
public double abs() { | |
return Math.hypot(re, im); | |
} | |
/** | |
* Returns the phase of this complex number. | |
* This quantity is also known as the <em>angle</em> or <em>argument</em>. | |
* | |
* @return the phase of this complex number, a real number between -pi and pi | |
*/ | |
public double phase() { | |
return Math.atan2(im, re); | |
} | |
/** | |
* Returns the sum of this complex number and the specified complex number. | |
* | |
* @param that the other complex number | |
* @return the complex number whose value is {@code (this + that)} | |
*/ | |
public Complex plus(Complex that) { | |
double real = this.re + that.re; | |
double imag = this.im + that.im; | |
return new Complex(real, imag); | |
} | |
/** | |
* Returns the result of subtracting the specified complex number from | |
* this complex number. | |
* | |
* @param that the other complex number | |
* @return the complex number whose value is {@code (this - that)} | |
*/ | |
public Complex minus(Complex that) { | |
double real = this.re - that.re; | |
double imag = this.im - that.im; | |
return new Complex(real, imag); | |
} | |
/** | |
* Returns the product of this complex number and the specified complex number. | |
* | |
* @param that the other complex number | |
* @return the complex number whose value is {@code (this * that)} | |
*/ | |
public Complex times(Complex that) { | |
double real = this.re * that.re - this.im * that.im; | |
double imag = this.re * that.im + this.im * that.re; | |
return new Complex(real, imag); | |
} | |
/** | |
* Returns the product of this complex number and the specified scalar. | |
* | |
* @param alpha the scalar | |
* @return the complex number whose value is {@code (alpha * this)} | |
*/ | |
public Complex scale(double alpha) { | |
return new Complex(alpha * re, alpha * im); | |
} | |
/** | |
* Returns the product of this complex number and the specified scalar. | |
* | |
* @param alpha the scalar | |
* @return the complex number whose value is {@code (alpha * this)} | |
* @deprecated Replaced by {@link #scale(double)}. | |
*/ | |
@Deprecated | |
public Complex times(double alpha) { | |
return new Complex(alpha * re, alpha * im); | |
} | |
/** | |
* Returns the complex conjugate of this complex number. | |
* | |
* @return the complex conjugate of this complex number | |
*/ | |
public Complex conjugate() { | |
return new Complex(re, -im); | |
} | |
/** | |
* Returns the reciprocal of this complex number. | |
* | |
* @return the complex number whose value is {@code (1 / this)} | |
*/ | |
public Complex reciprocal() { | |
double scale = re*re + im*im; | |
return new Complex(re / scale, -im / scale); | |
} | |
/** | |
* Returns the real part of this complex number. | |
* | |
* @return the real part of this complex number | |
*/ | |
public double re() { | |
return re; | |
} | |
/** | |
* Returns the imaginary part of this complex number. | |
* | |
* @return the imaginary part of this complex number | |
*/ | |
public double im() { | |
return im; | |
} | |
/** | |
* Returns the result of dividing the specified complex number into | |
* this complex number. | |
* | |
* @param that the other complex number | |
* @return the complex number whose value is {@code (this / that)} | |
*/ | |
public Complex divides(Complex that) { | |
return this.times(that.reciprocal()); | |
} | |
/** | |
* Returns the complex exponential of this complex number. | |
* | |
* @return the complex exponential of this complex number | |
*/ | |
public Complex exp() { | |
return new Complex(Math.exp(re) * Math.cos(im), Math.exp(re) * Math.sin(im)); | |
} | |
/** | |
* Returns the complex sine of this complex number. | |
* | |
* @return the complex sine of this complex number | |
*/ | |
public Complex sin() { | |
return new Complex(Math.sin(re) * Math.cosh(im), Math.cos(re) * Math.sinh(im)); | |
} | |
/** | |
* Returns the complex cosine of this complex number. | |
* | |
* @return the complex cosine of this complex number | |
*/ | |
public Complex cos() { | |
return new Complex(Math.cos(re) * Math.cosh(im), -Math.sin(re) * Math.sinh(im)); | |
} | |
/** | |
* Returns the complex tangent of this complex number. | |
* | |
* @return the complex tangent of this complex number | |
*/ | |
public Complex tan() { | |
return sin().divides(cos()); | |
} | |
/** | |
* Unit tests the {@code Complex} data type. | |
* | |
* @param args the command-line arguments | |
*/ | |
public static void main(String[] args) { | |
Complex a = new Complex(5.0, 6.0); | |
Complex b = new Complex(-3.0, 4.0); | |
System.out.println("a = " + a); | |
System.out.println("b = " + b); | |
System.out.println("Re(a) = " + a.re()); | |
System.out.println("Im(a) = " + a.im()); | |
System.out.println("b + a = " + b.plus(a)); | |
System.out.println("a - b = " + a.minus(b)); | |
System.out.println("a * b = " + a.times(b)); | |
System.out.println("b * a = " + b.times(a)); | |
System.out.println("a / b = " + a.divides(b)); | |
System.out.println("(a / b) * b = " + a.divides(b).times(b)); | |
System.out.println("conj(a) = " + a.conjugate()); | |
System.out.println("|a| = " + a.abs()); | |
System.out.println("tan(a) = " + a.tan()); | |
} | |
} | |
/****************************************************************************** | |
* Copyright 2002-2016, Robert Sedgewick and Kevin Wayne. | |
* | |
* This file is part of algs4.jar, which accompanies the textbook | |
* | |
* Algorithms, 4th edition by Robert Sedgewick and Kevin Wayne, | |
* Addison-Wesley Professional, 2011, ISBN 0-321-57351-X. | |
* http://algs4.cs.princeton.edu | |
* | |
* | |
* algs4.jar is free software: you can redistribute it and/or modify | |
* it under the terms of the GNU General Public License as published by | |
* the Free Software Foundation, either version 3 of the License, or | |
* (at your option) any later version. | |
* | |
* algs4.jar is distributed in the hope that it will be useful, | |
* but WITHOUT ANY WARRANTY; without even the implied warranty of | |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
* GNU General Public License for more details. | |
* | |
* You should have received a copy of the GNU General Public License | |
* along with algs4.jar. If not, see http://www.gnu.org/licenses. | |
******************************************************************************/ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
public class CrossCorrelation { | |
public static void main(String[] args) { | |
double[] source = {1, 2, 3, 4, 5, 6, 7, 8}; | |
double[] target = {1, 2, 3, 0, 0, 0, 0, 0}; | |
int n = source.length; | |
Complex[] sourceComplex = new Complex[n]; | |
Complex[] targetComplex = new Complex[n]; | |
for (int i = 0; i < n; i++) { | |
sourceComplex[i] = new Complex(source[i], 0); | |
targetComplex[i] = new Complex(target[i], 0); | |
} | |
Complex[] fftS = FFT.fft(sourceComplex); | |
Complex[] fftT = FFT.fft(targetComplex); | |
for (int i = 0; i < fftS.length; i++) { | |
fftS[i] = fftS[i].conjugate(); | |
} | |
Complex[] timeProduct = new Complex[fftS.length]; | |
for (int i = 0; i < fftS.length; i++) { | |
timeProduct[i] = fftS[i].times(fftT[i]); | |
} | |
Complex[] y = FFT.ifft(timeProduct); | |
System.out.println("id:\txcorr abs:\t\t\txcorr complexes:"); | |
for (int i = 0; i < y.length; i++) { | |
Complex c = y[i]; | |
System.out.println(i + "\t" +c.abs() + "\t\t\t" + c.toString()); | |
} | |
System.out.println("Max arg: " + argmax(y)); //I suppose it should equals 0 | |
} | |
public static int argmax(Complex[] a) | |
{ | |
double y = Double.MIN_VALUE; | |
int idx = -1; | |
for(int x = 0; x < a.length; x++) | |
{ | |
if(a[x].abs() > y) | |
{ | |
y = a[x].abs(); | |
idx = x; | |
} | |
} | |
return idx; | |
} | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/****************************************************************************** | |
* Compilation: javac FFT.java | |
* Execution: java FFT n | |
* Dependencies: Complex.java | |
* | |
* Compute the FFT and inverse FFT of a length n complex sequence. | |
* Bare bones implementation that runs in O(n log n) time. Our goal | |
* is to optimize the clarity of the code, rather than performance. | |
* | |
* Limitations | |
* ----------- | |
* - assumes n is a power of 2 | |
* | |
* - not the most memory efficient algorithm (because it uses | |
* an object type for representing complex numbers and because | |
* it re-allocates memory for the subarray, instead of doing | |
* in-place or reusing a single temporary array) | |
* | |
* | |
* % java FFT 4 | |
* x | |
* ------------------- | |
* -0.03480425839330703 | |
* 0.07910192950176387 | |
* 0.7233322451735928 | |
* 0.1659819820667019 | |
* | |
* y = fft(x) | |
* ------------------- | |
* 0.9336118983487516 | |
* -0.7581365035668999 + 0.08688005256493803i | |
* 0.44344407521182005 | |
* -0.7581365035668999 - 0.08688005256493803i | |
* | |
* z = ifft(y) | |
* ------------------- | |
* -0.03480425839330703 | |
* 0.07910192950176387 + 2.6599344570851287E-18i | |
* 0.7233322451735928 | |
* 0.1659819820667019 - 2.6599344570851287E-18i | |
* | |
* c = cconvolve(x, x) | |
* ------------------- | |
* 0.5506798633981853 | |
* 0.23461407150576394 - 4.033186818023279E-18i | |
* -0.016542951108772352 | |
* 0.10288019294318276 + 4.033186818023279E-18i | |
* | |
* d = convolve(x, x) | |
* ------------------- | |
* 0.001211336402308083 - 3.122502256758253E-17i | |
* -0.005506167987577068 - 5.058885073636224E-17i | |
* -0.044092969479563274 + 2.1934338938072244E-18i | |
* 0.10288019294318276 - 3.6147323062478115E-17i | |
* 0.5494685269958772 + 3.122502256758253E-17i | |
* 0.240120239493341 + 4.655566391833896E-17i | |
* 0.02755001837079092 - 2.1934338938072244E-18i | |
* 4.01805098805014E-17i | |
* | |
******************************************************************************/ | |
/** | |
* The {@code FFT} class provides methods for computing the | |
* FFT (Fast-Fourier Transform), inverse FFT, linear convolution, | |
* and circular convolution of a complex array. | |
* <p> | |
* It is a bare-bones implementation that runs in <em>n</em> log <em>n</em> time, | |
* where <em>n</em> is the length of the complex array. For simplicity, | |
* <em>n</em> must be a power of 2. | |
* Our goal is to optimize the clarity of the code, rather than performance. | |
* It is not the most memory efficient implementation because it uses | |
* objects to represents complex numbers and it it re-allocates memory | |
* for the subarray, instead of doing in-place or reusing a single temporary array. | |
* | |
* <p> | |
* For additional documentation, see <a href="http://algs4.cs.princeton.edu/99scientific">Section 9.9</a> of | |
* <i>Algorithms, 4th Edition</i> by Robert Sedgewick and Kevin Wayne. | |
* | |
* @author Robert Sedgewick | |
* @author Kevin Wayne | |
*/ | |
public class FFT { | |
private static final Complex ZERO = new Complex(0, 0); | |
// Do not instantiate. | |
private FFT() { } | |
/** | |
* Returns the FFT of the specified complex array. | |
* | |
* @param x the complex array | |
* @return the FFT of the complex array {@code x} | |
* @throws IllegalArgumentException if the length of {@code x} is not a power of 2 | |
*/ | |
public static Complex[] fft(Complex[] x) { | |
int n = x.length; | |
// base case | |
if (n == 1) { | |
return new Complex[] { x[0] }; | |
} | |
// radix 2 Cooley-Tukey FFT | |
if (n % 2 != 0) { | |
throw new IllegalArgumentException("n is not a power of 2"); | |
} | |
// fft of even terms | |
Complex[] even = new Complex[n/2]; | |
for (int k = 0; k < n/2; k++) { | |
even[k] = x[2*k]; | |
} | |
Complex[] q = fft(even); | |
// fft of odd terms | |
Complex[] odd = even; // reuse the array | |
for (int k = 0; k < n/2; k++) { | |
odd[k] = x[2*k + 1]; | |
} | |
Complex[] r = fft(odd); | |
// combine | |
Complex[] y = new Complex[n]; | |
for (int k = 0; k < n/2; k++) { | |
double kth = -2 * k * Math.PI / n; | |
Complex wk = new Complex(Math.cos(kth), Math.sin(kth)); | |
y[k] = q[k].plus(wk.times(r[k])); | |
y[k + n/2] = q[k].minus(wk.times(r[k])); | |
} | |
return y; | |
} | |
/** | |
* Returns the inverse FFT of the specified complex array. | |
* | |
* @param x the complex array | |
* @return the inverse FFT of the complex array {@code x} | |
* @throws IllegalArgumentException if the length of {@code x} is not a power of 2 | |
*/ | |
public static Complex[] ifft(Complex[] x) { | |
int n = x.length; | |
Complex[] y = new Complex[n]; | |
// take conjugate | |
for (int i = 0; i < n; i++) { | |
y[i] = x[i].conjugate(); | |
} | |
// compute forward FFT | |
y = fft(y); | |
// take conjugate again | |
for (int i = 0; i < n; i++) { | |
y[i] = y[i].conjugate(); | |
} | |
// divide by n | |
for (int i = 0; i < n; i++) { | |
y[i] = y[i].scale(1.0 / n); | |
} | |
return y; | |
} | |
/** | |
* Returns the circular convolution of the two specified complex arrays. | |
* | |
* @param x one complex array | |
* @param y the other complex array | |
* @return the circular convolution of {@code x} and {@code y} | |
* @throws IllegalArgumentException if the length of {@code x} does not equal | |
* the length of {@code y} or if the length is not a power of 2 | |
*/ | |
public static Complex[] cconvolve(Complex[] x, Complex[] y) { | |
// should probably pad x and y with 0s so that they have same length | |
// and are powers of 2 | |
if (x.length != y.length) { | |
throw new IllegalArgumentException("Dimensions don't agree"); | |
} | |
int n = x.length; | |
// compute FFT of each sequence | |
Complex[] a = fft(x); | |
Complex[] b = fft(y); | |
// point-wise multiply | |
Complex[] c = new Complex[n]; | |
for (int i = 0; i < n; i++) { | |
c[i] = a[i].times(b[i]); | |
} | |
// compute inverse FFT | |
return ifft(c); | |
} | |
/** | |
* Returns the linear convolution of the two specified complex arrays. | |
* | |
* @param x one complex array | |
* @param y the other complex array | |
* @return the linear convolution of {@code x} and {@code y} | |
* @throws IllegalArgumentException if the length of {@code x} does not equal | |
* the length of {@code y} or if the length is not a power of 2 | |
*/ | |
public static Complex[] convolve(Complex[] x, Complex[] y) { | |
Complex[] a = new Complex[2*x.length]; | |
for (int i = 0; i < x.length; i++) | |
a[i] = x[i]; | |
for (int i = x.length; i < 2*x.length; i++) | |
a[i] = ZERO; | |
Complex[] b = new Complex[2*y.length]; | |
for (int i = 0; i < y.length; i++) | |
b[i] = y[i]; | |
for (int i = y.length; i < 2*y.length; i++) | |
b[i] = ZERO; | |
return cconvolve(a, b); | |
} | |
// display an array of Complex numbers to standard output | |
private static void show(Complex[] x, String title) { | |
StdOut.println(title); | |
StdOut.println("-------------------"); | |
for (int i = 0; i < x.length; i++) { | |
StdOut.println(x[i]); | |
} | |
StdOut.println(); | |
} | |
/*************************************************************************** | |
* Test client. | |
***************************************************************************/ | |
/** | |
* Unit tests the {@code FFT} class. | |
* | |
* @param args the command-line arguments | |
*/ | |
public static void main(String[] args) { | |
int n = Integer.parseInt(args[0]); | |
Complex[] x = new Complex[n]; | |
// original data | |
for (int i = 0; i < n; i++) { | |
x[i] = new Complex(i, 0); | |
x[i] = new Complex(StdRandom.uniform(-1.0, 1.0), 0); | |
} | |
show(x, "x"); | |
// FFT of original data | |
Complex[] y = fft(x); | |
show(y, "y = fft(x)"); | |
// take inverse FFT | |
Complex[] z = ifft(y); | |
show(z, "z = ifft(y)"); | |
// circular convolution of x with itself | |
Complex[] c = cconvolve(x, x); | |
show(c, "c = cconvolve(x, x)"); | |
// linear convolution of x with itself | |
Complex[] d = convolve(x, x); | |
show(d, "d = convolve(x, x)"); | |
} | |
} | |
/****************************************************************************** | |
* Copyright 2002-2016, Robert Sedgewick and Kevin Wayne. | |
* | |
* This file is part of algs4.jar, which accompanies the textbook | |
* | |
* Algorithms, 4th edition by Robert Sedgewick and Kevin Wayne, | |
* Addison-Wesley Professional, 2011, ISBN 0-321-57351-X. | |
* http://algs4.cs.princeton.edu | |
* | |
* | |
* algs4.jar is free software: you can redistribute it and/or modify | |
* it under the terms of the GNU General Public License as published by | |
* the Free Software Foundation, either version 3 of the License, or | |
* (at your option) any later version. | |
* | |
* algs4.jar is distributed in the hope that it will be useful, | |
* but WITHOUT ANY WARRANTY; without even the implied warranty of | |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
* GNU General Public License for more details. | |
* | |
* You should have received a copy of the GNU General Public License | |
* along with algs4.jar. If not, see http://www.gnu.org/licenses. | |
******************************************************************************/ |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment