Как извлечь пиковые значения из автокоррелированных данных в MATLAB?

Question

У меня есть информация (20 000 кадров данных) об аудиодорожке, которую я автоматически коррелировал с помощью:

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

И это выглядит так:

альтернативный текст http://a.imageshack.us/img809/3775/plot.jpg

Что, надеюсь, пока хорошо. В идеале я хотел бы иметь возможность взять номер кадра, который соответствует самой высокой части второго пика. Я прочитал и попробовал множество различных методов, но я просто не могу получить его, чтобы получить информацию для меня.

Кто-нибудь сможет пролить свет на то, что я должен делать?

Большое спасибо!

---

edit1: Я пытался использовать findpeaks, но мне это не помогло. Я не уверен, что это потому, что я использую неправильные данные или нет.

edit2: В настоящее время я тестирую метод для использования только с этой звуковой дорожкой, но вскоре я хочу расширить его, чтобы я мог выполнить этот метод для целого каталога файлов, поэтому я Нужен скрипт, который может обнаруживать пики, а не находить информацию сам.

edit3: Мой .M файл:

[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

РЕДАКТИРОВАТЬ :

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

и значения запаздывания, соответствующие отмеченным пикам:

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

Обратите внимание, чтомы сгладили сигнал до вычисления производной функции автокорреляции, чтобы найти точек экстремумов

Answer 2

Вы можете использовать findpeaks дважды.Во-первых, чтобы получить начальную оценку пиков и использовать полученные результаты для точной настройки входных параметров второго вызова findpeaks.Из вопроса кажется, что вы хотите рассчитать значение шага.Вот код:

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

Выше приведена модифицированная версия кода, заданного здесь .

Answer 3

Если у вас есть набор инструментов для обработки сигналов, я думаю, что функция findpeaks должна выполнить эту работу за вас.

Что-то вроде:

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

Я привык использовать его для данных интенсивности изображения, которые не имеют такой большой локальной дисперсии, какваши данные (Эти «толстые» участки вашего графика - это просто множество точек данных, быстро и быстро поднимающихся и опускающихся, верно?).Возможно, вам придется сначала сгладить данные, чтобы они заработали.

Answer 4

В качестве первого шага вы должны использовать второй выходной аргумент xcorr, чтобы получить правильные единицы измерения на графике:

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

Теперь вы можете использовать свой любимый метод, чтобы получить лаг для второго пика; кнопка проводника данных на рисунке подойдет, если вам просто нужно сделать это один раз. Если вам нужно сделать это несколько раз, я не вижу способа избежать попадания в findpeaks или подобное.

Answer 5

1. Сглаживание данных с использованием фильтра нижних частот (или усреднение каждой выборки по количеству окружающих выборок).
2. Найдите пик в центре, отыскивая наибольшее значение выборки.
3. найдите долину справа от пика путем поиска первого образца, который имеет более высокое значение, чем его предшественник.
4. Найдите пик справа от долины путем поиска первого образца в меньшем размерезначение, чем его предшественник.