I am trying to implement counting cross correlation in java based on its properties (the fifth one with FFT transform). It should work properly with signals of different length. I am using simple FFT library from Princeton University (FFT.java and Complex.java). Independently of the target signal I receive a cross correlation peak at 3rd index. For the example below it should be placed at 0 index because both signals are the same at first three places. I make that solution to have a possibility to synchronize two different size audio files (the one is a part of second file).
@Test
public void find_synchro_by_cross_correlation() {
double[] source = {1, 2, 3, 4, 5, 6, 7, 8};
double[] target = {1, 2, 3, 0, 0, 0, 0, 0}; //signal padded 0's
int n = source.length;
//creating complex arrays
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);
}
//counting FFT of both signals
Complex[] fftS = FFT.fft(sourceComplex);
Complex[] fftT = FFT.fft(targetComplex);
//conjugate the first one
for (int i = 0; i < fftS.length; i++) {
fftS[i] = fftS[i].conjugate();
}
// point-wise multiply
Complex[] timeProduct = new Complex[fftS.length];
for (int i = 0; i < fftS.length; i++) {
timeProduct[i] = fftS[i].times(fftT[i]);
}
// compute inverse FFT
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));
Assert.assertEquals(0, DSP.argmax(y)); //I assume that the peak should be at 0 index
}
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;
}
Invocation result:
id: xcorr abs: xcorr complexes:
0 14.0 14.0 + 5.551115123125783E-16i
1 16.0 16.0 + 1.133107779529596E-15i
2 26.0 26.0 - 8.881784197001252E-16i
3 44.0 44.0 - 2.021286199229721E-15i
4 38.0 38.0 - 5.551115123125783E-16i
5 32.0 32.0 - 6.432490598706545E-16i
6 26.0 26.0 + 8.881784197001252E-16i
7 20.0 20.0 + 1.5314274795707798E-15i
Max arg: 3 //I thought it should equals 0
For the refererence/debugging purposes here is a gist with sources to download full code.
That matlab/octave code makes the proper job (but still it doesn't use the FFT algorithms):
pkg load signal
x = [2 3 4];
y = 1:1:10;
nx = length(x);
ny = length(y);
cc = nan(1,ny-nx+1);
for ii = 0 : ny-nx
id = (1:nx) + ii;
cc(ii+1) = sum(x.*y(id))/(sqrt(sum(x.^2)*sum(y(id).^2)));
end
[ccmx,idmx] = max(cc)
There is another solution for that problem but still it takes so much time, complexity of the solution is n^2 lgn
Y = [2 3 4];
X = 1:1:10;
lngX = length(X);
lngY = length(Y);
assert(lngX >= lngY);
lags = 0:(lngX-lngY);
for i = lags
c(i+1) = xcorr(X(i+1:i+lngY) - mean(X(i+1:i+lngY)), Y - mean(Y),0,'coeff');
end
[m,i]=max(c);
printf('max=%f, lag=%d\n',c(i),lags(i));
plot(lags,c,'-',lags(i),c(i),'*r');
UPDATE
I've found kind of solution, but I am not satisfied its complexity. I rewrite last matlab example. I am counting FFT for the shorter signal only once as well as part of coefficiency (because of efficiency). I am counting the correlation for the 0 offset (the first element in coefficiency array) for every possible lag. The peak of this calculation shows how much signal has been lagged. If you would like to count FFT on arrays with length different that n=power2 then you have to change FFT library. You can use FFTW java wrapper or JTransforms. Below you can read code.
@Test
public void find_synchro_by_circular_convolution() {
double[] source = {1, 2, 3, 4, 5, 6, 7, 8};
double[] target = {4, 5, 6, 7};
int n = target.length;
int lags = source.length - n + 1;
//normalize(target);
Complex[] targetComplex = convertToComplex(target);
double sumsqTarget = sumsq(target);
Complex[] fftTarget = FFT.fft(targetComplex);
Complex[] c = new Complex[lags];
for (int i = 0; i < lags; i++) {
//double[] shrinkedSource = normalize(source, i, n);
double[] shrinkedSource = new double[n];
shrinkArray(source, shrinkedSource, i, n);
Complex[] sourceComplex = convertToComplex(shrinkedSource);
//cross correlation coefficient
Complex[] fftSource = FFT.fft(sourceComplex);
conjugate(fftSource);
multiplyComplexXInPlace(fftSource, fftTarget);
Complex[] y = FFT.ifft(fftSource);
double rms = rms(sumsqTarget, shrinkedSource);
c[i] = y[0].scale(rms);
}
System.out.println("id:\txcorr abs:\t\t\txcorr complexes:");
for (int i = 0; i < lags; i++) {
Complex ci = c[i];
System.out.println(i + "\t" +ci.abs() + "\t\t\t" + ci.toString());
}
System.out.println("Max arg: " + DSP.argmax(c));
}
@NotNull
private Complex[] convertToComplex(double[] target) {
int n = target.length;
Complex[] targetComplex = new Complex[n];
for (int i = 0; i < n; i++) {
targetComplex[i] = new Complex(target[i], 0);
}
return targetComplex;
}
private void multiplyComplexXInPlace(Complex[] x, Complex[] y) {
for (int j = 0; j < x.length; j++) {
x[j] = x[j].times(y[j]);
}
}
private void conjugate(Complex[] fftSource) {
for (int j = 0; j < fftSource.length; j++) {
fftSource[j] = fftSource[j].conjugate();
}
}
private double[] normalize(double[] source, int i, int n) {
double[] normalized = new double[n];
shrinkArray(source, normalized, i, n);
normalize(normalized);
return normalized;
}
private void normalize(double[] target) {
double meanTarget = mean(target);
for (int i = 0; i < target.length; i++) {
target[i] = target[i] - meanTarget;
}
}
private void shrinkArray(double[] source, double[] destination, int startPosition, int length) {
System.arraycopy(source, startPosition, destination, 0, length);
}
private double mean(double[] x) {
double sum = 0;
int n = x.length;
if(n != 0) {
for (int i = 0; i < n; i++) {
sum += x[i];
}
return sum / x.length;
}
return sum;
}
private double rms(double cached_sumsq_x, float[] y) {
double sumsq_x = cached_sumsq_x;
double sumsq_y = sumsq(y);
return 1/Math.sqrt(sumsq_x*sumsq_y);
}
private double rms(double cached_sumsq_x, double[] y) {
double sumsq_x = cached_sumsq_x;
double sumsq_y = sumsq(y);
return 1/Math.sqrt(sumsq_x*sumsq_y);
}
private double rms(float[] x, float[] y) {
double sumsq_x = sumsq(x);
double sumsq_y = sumsq(y);
return 1/Math.sqrt(sumsq_x*sumsq_y);
}
private double sumsq(float[] x) {
double sumsq = 0.0;
int n = x.length;
for (int i = 0; i < n; i++) {
sumsq += x[i]*x[i];
}
return sumsq;
}
private double sumsq(double[] x) {
double sumsq = 0.0;
int n = x.length;
for (int i = 0; i < n; i++) {
sumsq += x[i]*x[i];
}
return sumsq;
}
The result of code:
id: xcorr abs: xcorr complexes:
0 0.9759000729485332 0.9759000729485332
1 0.9941037221861875 0.9941037221861875
2 0.9990767240206739 0.9990767240206739
3 1.0 1.0
4 0.9995437747834998 0.9995437747834998
Max arg: 3
But... come on. Simple linear sollution is faster! Is it really so hard count it in frequency domain?
import static java.lang.Math.abs;
public class LinearSynchronisation {
public long findSynchroSamples(Float[] longFile, Float[] shortfile) {
double lowestDifference = Double.MAX_VALUE;
long maxCorrelationSampleNumer = 0;
for (int j = 0; j < longFile.length - shortfile.length; j++) {
double diff = getCorrelation(longFile, shortfile, j);
if (diff < lowestDifference) {
maxCorrelationSampleNumer = j;
lowestDifference = diff;
}
}
return maxCorrelationSampleNumer;
}
private double getCorrelation(Float[] longFile, Float[] shortFile, int start) {
double diff = 0;
for (int i = 0; i < shortFile.length; i++) {
diff += abs(longFile[start+i] - shortFile[i]);
}
return diff;
}
public static void main(String[] args) {
Float[] source = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f};
Float[] target = {4.0f, 5.0f, 6.0f, 7.0f};
LinearSynchronisation ls = new LinearSynchronisation();
long synchroSamples = ls.findSynchroSamples(source, target);
System.out.println(synchroSamples);
}
}
