Take the 2-minute tour ×
Signal Processing Stack Exchange is a question and answer site for practitioners of the art and science of signal, image and video processing. It's 100% free, no registration required.

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.

share|improve this question
add comment

migrated from electronics.stackexchange.com Nov 30 '11 at 18:18

This question came from our site for electronics and electrical engineering professionals, students, and enthusiasts.

4 Answers

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 (element-wise multiplication)
  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:

peak in 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 -1 and 1 (a negative value indicating one of the signals decreases as the other increases) 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.

share|improve this answer
1  
The zero padding should be at least N = size(a)+size(b)-1, preferably rounded up to a power of 2. To get a value between -1 and 1, divide by norm(a)*norm(b), which gives the cosine of the angle between the two vectors in N-space for the given lag (i.e. circular shift modulo N). At the extreme lags, there aren't many overlapping samples (only one at the far extreme), so dividing by norm(a)*norm(b) will bias these correlations toward 0 (i.e. showing their relative orthogonality in N-space). –  eryksun Oct 30 '10 at 1:41
1  
I think there may be an error in the description. Shouldn't multiplying the FFTs together term by term give the FFT of the convolution of the signals, not the FFT of the cross-correlation? As I understand it, to get the FFT of the cross-correlation, it is necessary to use the complex conjugate of one of the FFT vectors in the term-by-term multiplications before taking the iFFT. –  Dilip Sarwate Nov 30 '11 at 20:36
    
@DilipSarwate: Yes, you're right. You can also reverse one signal in the time direction, which I added to the answer. –  endolith Nov 30 '11 at 20:49
1  
"Why is time reversal hard to do in hardware?" In lots of cases, data are stored in systolic arrays in the expectation that computations are local, i.e. $x[i]$, stored in the $i$-th cell, interacts only with its nearest neighbors $x[\pm i]$. Sending $x[i]$ to cell #$(N-i)$ and sending $x[N-i]$ to cell #$i$, and doing this for all $i$ increasing wiring costs, wiring delays (and hence reduces maximum achievable clock rate), and also, because all the wires must cross over each other, creates routing problems. It should be avoided if possible, and in this case, it is avoidable. –  Dilip Sarwate Nov 30 '11 at 22:05
1  
@Leo element-wise multiplication. n-by-1 array x n-by-1 array = n-by-1 array I called this "sample-by-sample" in the answer. –  endolith Jan 22 at 20:34
show 6 more comments

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.

share|improve this answer
1  
You mentioned that 0 is no correlation and 1 is total correlations. I just want to note that -1 is complete negatively correlated. As in, -1 implies that sample 1 is the opposite of sample 2. If you think about it on an X,Y graph, it is a line with positive slope versus a line with negative slope. And as you get closer to 0 the line gets "fatter". –  Kellenjb Oct 29 '10 at 13:23
    
@kellenjb, Yes, but I would probably word it, the magnitude of correlation what you are probably interested in. a 1 or a -1 mean the signals directly affect each other. –  Kortuk Oct 29 '10 at 14:52
add comment

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))
share|improve this answer
    
Does this mean the final size of the array should be 2*size (a) + size(b) - 1 or 2*size (b) + size (a) - 1? But in either case the two padded arrays are of different sizes. What is the consequence of padding with too many zeros? –  RobertKJ Nov 29 '11 at 1:20
    
@RobertK The cross-correlation array needs to be of length at least the sum of lengths of a and b (minus one) as eryksun says in his answer. For simplicity, the length is often taken to be twice the length of the longer vector (sometimes rounded up to the next larger power of $2$ in order to use an efficient FFT). The choice helps when the customer belatedly decides he also wants the autocorrelation of the longer vector. One consequence of padding with too many zeroes is additional computation but this might be ameliorated by more efficient FFT implementations. –  Dilip Sarwate Nov 30 '11 at 20:54
    
@RobertKJ: You're sliding b along a, with one output per shift, a minimum overlap of one sample. That yields size(a) positive lags and size(b) - 1 negative lags. Using the inverse transform of the product of N-point DFTs, indices 0 through size(a)-1are the positive lags, and indices N-size(b)+1 through N-1 are the negative lags in reverse order. –  eryksun Nov 30 '11 at 23:56
add comment

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.

correlation

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.