Geek Answers Handbook

How can I extract peak values ​​from autocorrelated data in MATLAB?

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

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

And it looks like this:

alt text http://a.imageshack.us/img809/3775/plot.jpg

Which, I hope, is so good so far. Ideally, I would like to have a frame number that corresponds to the highest part of the second peak. I read and tried to load different methods, but I just can't get it to get the information for me.

Can anyone shed some light on what I have to do?

Many thanks!

edit1: I tried using findpeaks , but it doesn't seem to work for me. I am not sure if this is because I am using the wrong data or not.

edit2: I am currently testing a method that can only be used on this sound track, but I want to expand it soon so that I can execute this method in the whole file directory, so I am kind of needing a script that can detect peaks rather than find information itself.

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
+4
source share
5 answers

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')

alt text

and lag values ​​corresponding to the marked peaks:

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

Note that we smoothed the signal before calculating the derivative of the autocorrelation function to find extreme points

+5
source
  • Smooth data using a low-pass filter (or average each sample with a number of surrounding samples).
  • Find the peak in the center by looking for the highest sample value.
  • find the valley to the right of the peak, looking for the first sample that has a higher value than its predecessor.
  • Find the peak to the right of the valley, looking for the first sample at a lower value than its predecessor.
+1
source

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

Sort of:

 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 use it for image intensity data that does not have the same local variance as your data (are these “thick” sectors of your plot just a few data points that grow quickly and quickly?). You may need to smooth the data a bit to make it work.

+1
source

As a first step, you should use the second output argument xcorr to properly configure your units on the chart:

 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 delay 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 this several times, I see no way to avoid getting into findpeaks or the like.

+1
source

You can use findpeaks twice. First, to get an initial peak estimate and use this 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 version is a modified version of the code here .

+1
source

More articles:


All Articles