How can I extract peak values from auto-correlated data in MATLAB?

Question

I have information (20,000 frames of data) about an audio track that I have auto-correlated using:

[r,lags] = xcorr(XX,XX,'biased');

And it looks like this:

Which hopefully is so far so good. Ideally I would like to be able to take the frame number that corresponds to the highest part of the second peak. I've read around and tried a load of different methods but I just can't seem to get it to retrieve the information for me.

Would anybody be able to shed some light on what I have to do?

Many thanks!

---

edit1: I have tried using findpeaks, but it doesn't seem to work for me. I'm not sure if that's because I'm using the wrong data or not.

edit2: I'm currently testing a method to use on just this audio track, but soon I want to expand it so that I can perform this method on a whole directory of files, so I kind of need a script that can detect peaks rather than finding the information myself.

edit3: My .M file:

[y, fs, nb] = wavread('Three.wav');                 %# Load the signal into variable y

frameWidth = 441;                                   %# 10ms
numSamples = length(y);                             %# Number of samples in y
numFrames = floor(numSamples/frameWidth);           %# Number of full frames in y
energy = zeros(1,numFrames);                        %# Initialize energy
startSample = zeros(1,numFrames);                   %# Initialize start indices
endSample = zeros(1,numFrames);                     %# Initialize end indices

for frame = 1:numFrames                             %# Loop over frames
  startSample(frame) = (frame-1)*frameWidth+1;      %# Starting index of frame
  endSample(frame) = frame*frameWidth;              %# Ending index of frame
  frameIndex = startSample(frame):endSample(frame); %# Indices of frame samples
  energy(frame) = sum(y(frameIndex).^2);            %# Calculate frame energy
end                                                 %# End loop

XX = filtfilt(ones(1,10)/10, 1, energy);            %# Smooths signal

[r,lags] = xcorr(XX,XX,'biased');                   %# Auto-correlates the data
plot(lags,r), xlabel('lag'), ylabel('xcorr')        %# Plots data

Answer 1

EDIT:

%# load the signal
[y, fs, nb] = wavread('Three.wav');
y = mean(y,2);                               %# stereo, take avrg of 2 channels

%# Calculate frame energy
fWidth = round(fs*10e-3);                    %# 10ms
numFrames = floor(length(y)/fWidth);
energy = zeros(1,numFrames);
for f=1:numFrames
  energy(f) = sum( y((f-1)*fWidth+1:f*fWidth).^2 );
end

%# smooth the signal (moving average with window size = 1% * length of data)
WINDOW_SIZE = round(length(energy) * 0.01);  %# 200
XX = filtfilt(ones(1,WINDOW_SIZE)/WINDOW_SIZE, 1, energy);

%# auto-correlation
[r,lags] = xcorr(XX, 'biased');

%# find extrema points
dr = diff(r);
eIdx = find(dr(1:end-1) .* dr(2:end) <= 0) + 1;

[~,loc] = sort(r(eIdx), 'descend');
loc = loc(1:min(3,end));                     %# take the highest 3 values

%# plot
plot(lags,r), hold on
plot(lags(eIdx), r(eIdx), 'g*')
plot(lags(eIdx(loc)), r(eIdx(loc)), 'ro')
hold off, xlabel('lag'), ylabel('xcorr')

and the lag values corresponding to the marked peaks:

>> lags( eIdx(loc) )
ans =
           0       -6316        6316

Note that we smoothed the signal prior to computing the derivative of the autocorrelation function in order to find the extrema points

Answer 2

1. Smooth the data using a lowpass filter (or average each sample with a number of surrounding samples).
2. Find the peak in the center by looking for the highest sample value.
3. find the valley to the right of the peak by searching for the first sample that has a higher value than its predecessor.
4. Find the peak to the right of the valley by searching for the first sample to a smaller value than its predecessor.

Answer 3

If you have the signal processing toolbox, I think the findpeaks function should do the job for you.

Something like:

th = x //(some value that will overlook noise in the data. See the documentation)
[peaks locs] = findpeaks(a,'threshold',th)

http://www.mathworks.com/access/helpdesk/help/toolbox/signal/findpeaks.html

I'm used to using it on image intensity data, which doesn't have quite as much local variance as your data (Those "thick" looking sections of your plot are just many data points going up and down quickly, right?). You might need to smooth the data out a bit first to get it to work.

Answer 4

As a first step, you should use the second output argument of xcorr to get your units on the plot correct:

Fs = length(data)./T; % or replace with whatever your sample frequency is (44.1 kHz?)
[a,lags] = xcorr(data,data,'biased');
plot(lags./Fs,a);

Now you can use your favorite method to get the lag for the second peak; the data explorer button in the figure will do if you just need to do it once. If you need to do it repeatedly, I don't see any way to avoid getting into findpeaks or similar.

Answer 5

You can use findpeaks twice. First to get an initial estimate of the peaks and use that results to fine tune the input parameters of the second call to findpeaks. From the question it seems that you want to calculate the pitch value. Here is the code :

maxlag = fs/50;
r = xcorr(x, maxlag, 'coeff');
r_slice = r(ceil(length(r)/2) : length(r));
[pksh,lcsh] = findpeaks(r_slice);

if length(lcsh) > 1
short = mean(diff(lcsh));
else
    short = lcsh(1)-1;
end

[pklg,lclg] = findpeaks(r_slice,'MinPeakDistance',ceil(short),'MinPeakheight',0.3);

if length(lclg) > 1
    long = mean(diff(lclg));
else
    if length(lclg) > 0
        long = lclg(1)-1;
    else
        long = -1;
    end
end

if long > 0
    pitch = fs / long; % fs is sample rate in Hertz

The above is a modified version of the code given here.