如何从 MATLAB 中的自相关数据中提取峰值？

Question

我有关于音轨的信息（20,000 帧数据），我使用它自动关联：

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

它看起来像这样：

希望到目前为止一切顺利。理想情况下，我希望能够获取与第二个峰值的最高部分相对应的帧号。我已经阅读并尝试了很多不同的方法，但我似乎无法让它为我检索信息。

有人能够阐明我必须做什么吗？

非常感谢！

---

编辑1： 我尝试过使用 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

1. 使用低通滤波器平滑数据（或对每个样本与多个周围样本进行平均）。
2. 通过寻找最高样本值来找到中心的峰值。
3. 通过搜索比其前一个值更高的第一个样本，找到峰值右侧的谷值。
4. 通过搜索比其前一个值更小的第一个样本，找到谷右侧的峰值。

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

您可以使用 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

上面是此处给出的代码的修改版本。