How do I implement cross-correlation to prove two audio files are similar?
Cross-correlation and convolution are closely related.  In short, to do convolution with FFTs, you 


zero-pad the input signals (add zeros to the end so that at least half of the wave is "blank")
take the FFT of both signals
multiply the results together (element-wise multiplication)
do the inverse FFT


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.