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.
Cross-correlation and convolution are closely related. In short, to do convolution with FFTs, you
a and b (add zeros to the end of each. The zero padding should fill the vectors until they reach a size of at least N = size(a)+size(b)-1)conv(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros))
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(a, b) = ifft(fft(a_and_zeros) * fft(b_and_zeros[reversed]))corr(a, b) = ifft(fft(a_and_zeros) * conj(fft(b_and_zeros)))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.
If the signals are real, you can use real FFTs (RFFT/IRFFT) and save half your computation time by only calculating half of the spectrum.
Also you can save computation time by padding to a larger size that the FFT is optimized for (such as a 5-smooth number for FFTPACK, a ~13-smooth number for FFTW, or a power of 2 for a simple hardware implementation).
Here's an example in Python of FFT correlation compared with brute-force correlation: https://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:
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.
mode='full' they should do the same thing
$\endgroup$
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.
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))
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?
$\endgroup$
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.
$\endgroup$
if you are using Matlab try the cross correlate function:
c= xcorr(x,y)
Here's the Matlab documentation:
xcorrestimates 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, wherexandyare lengthNvectors (N > 1). Ifxandyare not the same length, the shorter vector is zero-padded to the length of the longer vector.correlation http://www.mathworks.com/help/toolbox/signal/ref/eqn1263487323.gif
A quick and simple way to compare audio files. Take the audio file, make a copy, in a daw, paste them side by side, in 2 stereo channels, invert the phase on one of the stereo tracks, align both of the file at the beginning in zoom mode, be sure that the both files have same amplitude at the beginning, then play ,if there's total silence, then both files are identical, if there's a difference you'll hear it pretty clearly!.
The easiest way to find the difference, IMO, is to subtract the two audio signals in the time domain. If they are equal, the result at every time point will be zero. If they are not equal, the difference between them will be left after subtraction and you can listen to it directly. A quick measure of how similar they are would be the RMS value of this difference. This is often done in audio mixing and mastering to hear the difference of an MP3 vs WAV file for example. (Inverting the phase of one signal and adding them is the same as subtracting. This is the method used when this is done in DAW software.) They must be perfectly time aligned for this to work. If they are not you could develop an algorithm to align them such as detecting the top ten peaks, calculating the average offset of the peaks, and shifting one signal.
Transforming to the frequency domain and comparing the power spectra of the signals like you propose is ignoring some time domain information. For example audio played in reverse would have the same spectrum when played forward. Thus, two very different audio signals could have the exact same spectrum.
As most here wrote you should use correlation.
Just take 2 factors under consideration:
