How do I implement cross-correlation to prove two audio files are similar?

Question

I have to do cross correlation of two audio file to prove they are similar. I have taken the FFT of the two audio files and have their power spectrum values in separate arrays.

How should I proceed further to cross-correlate them and prove that they're similar? Is there a better way to do it? Any basic ideas will be helpful for me to learn and apply it.

Answer 1

Cross-correlation and convolution are closely-related. In short, to do convolution with FFTs, you

1. zero-pad the input signals (add zeros to the end so that at least half of the wave is "blank")
2. take the FFT of both signals
3. multiply the results together (sample-by-sample)
4. do the inverse FFT

conv = ifft(fft(a and zeroes) * fft(b and zeroes))

You need to do the zero-padding because the FFT method is actually circular cross-correlation, meaning the signal wraps around at the ends. So you add enough zeros to get rid of the overlap, to simulate a signal that is zero out to infinity.

To get cross-correlation instead of convolution, you either need to time-reverse one of the signals before doing the FFT, or take the complex conjugate of one of the signals after the FFT:

• corr = ifft(fft(a and zeroes) * fft(b and zeroes[reversed]))
• corr = ifft(fft(a and zeroes) * conj(fft(b and zeroes)))

whichever is easier with your hardware/software. For autocorrelation (cross-correlation of a signal with itself), it's better to do the complex conjugate, because then you only need to calculate the FFT once.

Here's an example in Python of FFT correlation compared with brute-force correlation: http://stackoverflow.com/a/1768140/125507

This will give you the cross-correlation function, which is a measure of similarity vs offset. To get the offset at which the waves are "lined up" with each other, there will be a peak in the correlation function:

http://jfk-records.com/NRC_Science/science.htm

The x value of the peak is the offset, which could be negative or positive.

I've only seen this used to find the offset between two waves. You can get a more precise estimate of the offset (better than the resolution of your samples) by using parabolic/quadratic interpolation on the peak.

To get a similarity value between 0 and 1 you'd need to scale the amplitude according to the length of the inputs, length of the FFT, your particular FFT implementation's scaling, etc. The autocorrelation of a wave with itself will give you the value of the maximum possible match.

Note that this will only work on waves that have the same shape. If they've been sampled on different hardware or have some noise added, but otherwise still have the same shape, this comparison will work, but if the wave shape has been changed by filtering or phase shifts, they may sound the same, but won't correlate as well.

Answer 2

Correlation is a way to express the similarity of two timeseries (audio samples in your case) in one number. It is an adaptation of covariance which is implemented as follows:

period = 1/sampleFrequency;
covariance=0;

for (iSample = 0; iSample<nSamples; iSample++)
    covariance += (timeSeries_1(iSample)*timeSeries_2(iSample))/period;
    //Dividing by `period` might not even be necessary

The correlation is the normalized version of covariance, which is the covariance divided by the product of the standard deviations of both the time series. The correlation will yield a 0 when there is no correlation (totally not similar) and a 1 for total correlation (totally similar).

You can imagine that two sound samples might be similar but are not synchronized. That's where cross correlation comes in. You calculate the correlation between the time series where you have one of them shifted by one sample:

for (iShift=0; iShift<nSamples; iShift++)
    xcorr(iShift) = corr(timeSeries_1, timeSeries_2_shifted_one_sample);

Then seek out the maximum value in the corr series and you're done. (or stop if you've found a sufficient correlation) There's slightly more to it of course. You must implement the standard deviation and you must do some memory management and implement the time shifting stuff. If all your audio samples are equal in length, you might do without normalizing the covariance and go ahead and calculate the cross-covariance.

A cool relation to your earlier question: Fourier analysis is just an adaptation of the cross covariance. Rather than shifting one time series and calculating the covariances with the other signal, you calculate the covariances between one signal and a number of (co)sine waves with different frequencies. It's all based on the same principle.

Answer 3

In signal processing the cross-correlation (xcorr in MATLAB) is a convolution operation with one of the two sequences reversed. Since time reversal corresponds to complex conjugation in the frequency domain, you can use the DFT to compute the cross-correlation as follows:

R_xy = ifft(fft(x,N) * conj(fft(y,N)))

where N = size(x) + size(y) - 1 (preferably rounded up to a power of 2) is the length of the DFT.

Multiplication of DFTs is equivalent to circular convolution in time. Zero padding both vectors to length N keeps the circularly shifted components of y from overlapping with x, which makes the result identical to the linear convolution of x and time reversed y.

A lag of 1 is a right circular shift of y, while a lag of -1 is a left circular shift. The cross-correlation is simply the sequence of dot products for all lags. Based on standard fft ordering, these will be in an array that can be accessed as follows. Indices 0 through size(x)-1 are the positive lags. Indices N-size(y)+1 to N-1 are the negative lags in reverse order. (In Python the negative lags can be accessed conveniently with negative indices such as R_xy[-1].)

You can think of the zero-padded x and y as N-dimensional vectors. The dot product of x and y for a given lag is |x|*|y|*cos(theta). The norms of x and y are constant for circular shifts, so dividing them out leaves just the varying cosine of the angle theta. If x and y (for a given lag) are orthogonal in N-space, the correlation is 0 (i.e. theta = 90 degrees). If they're co-linear, the value is either 1 (positively correlated) or -1 (negatively correlated, i.e. theta = 180 degrees). This leads to the cross-correlation normalized to unity:

R_xy = ifft(fft(x,N) * conj(fft(y,N))) / (norm(x) * norm(y))

This can be made unbiased by recomputing the norms for just the overlapping parts, but then you may as well do the entire computation in the time domain. Also, you'll see different versions of normalization. Instead of being normalized to unity, sometimes the cross-correlation is normalized by M (biased), where M = max(size(x), size(y)), or M-|m| (an unbiased estimate of the mth lag).

For maximum statistical significance the mean (DC bias) should be removed before computing the correlation. This is called the cross-covariance (xcov in MATLAB):

x2 = x - mean(x)
y2 = y - mean(y)
phi_xy = ifft(fft(x2,N) * conj(fft(y2,N))) / (norm(x2) * norm(y2))

Answer 4

if you are using Matlab try the cross correlate function:

c= xcorr(x,y)

Here's the Matlab documentation:

xcorr estimates the cross-correlation sequence of a random process. Autocorrelation is handled as a special case.

...

c = xcorr(x,y) returns the cross-correlation sequence in a length 2*N-1 vector, where x and y are length N vectors (N > 1). If x and y are not the same length, the shorter vector is zero-padded to the length of the longer vector.