嘈杂正弦时间序列中的实时峰值检测 [英] Real-time peak detection in noisy sinusoidal time-series

问题描述

我一直在尝试实时实时检测正弦时间序列数据中的峰值，但是到目前为止，我还没有成功.我似乎找不到一个实时算法来以合理的准确度检测正弦信号中的峰值.我要么没有检测到峰值，要么在正弦波上得到了数以十万计的点被检测为峰值.

I have been attempting to detect peaks in sinusoidal time-series data in real time, however I've had no success thus far. I cannot seem to find a real-time algorithm that works to detect peaks in sinusoidal signals with a reasonable level of accuracy. I either get no peaks detected, or I get a zillion points along the sine wave being detected as peaks.

对于类似于正弦波并可能包含一些随机噪声的输入信号，有什么好的实时算法?

作为一个简单的测试用例，考虑一个固定的正弦波，其频率和幅度始终相同. (确切的频率和幅度无关紧要；我以8 KS/s的采样率任意选择了60 Hz的频率，+/-1的幅度.)下面的MATLAB代码将生成这样的正弦信号:

As a simple test case, consider a stationary, sine wave that is always the same frequency and amplitude. (The exact frequency and amplitude don't matter; I have arbitrarily chosen a frequency of 60 Hz, an amplitude of +/− 1 unit, at a sampling rate of 8 KS/s.) The following MATLAB code will generate such a sinusoidal signal:

dt = 1/8000;
t  = (0:dt:(1-dt)/4)';
x  = sin(2*pi*60*t);

使用 Jean-Paul开发和发布的算法，我要么没有检测到峰(左)，要么有几千个峰值" (正确)检测到:

Using the algorithm developed and published by Jean-Paul, I either get no peaks detected (left) or a zillion "peaks" detected (right):

在经验法则"让·保罗(Jean-Paul)给出了，但是到目前为止，我一直无法获得预期的结果.

I've tried just about every combination of values for these 3 parameters that I could think of, following the "rules of thumb" that Jean-Paul gives, but I have so far been unable to get my expected result.

我发现了由Eli Billauer开发和发布的另一种算法 会给我想要的结果吗?例如:

I found an alternative algorithm, developed and published by Eli Billauer, that does give me the results that I want—e.g.:

即使Eli Billauer的算法简单得多，并且确实可以可靠地产生我想要的结果，但它不适用于实时应用程序.

Even though Eli Billauer's algorithm is much simpler and does tend to reliably produce the results that I want, it is not suitable for real-time applications.

作为我想应用这种算法的信号的另一个示例，请考虑Eli Billauer针对他自己的算法给出的测试用例:

As another example of a signal that I'd like to apply such an algorithm to, consider the test case given by Eli Billauer for his own algorithm:

t = 0:0.001:10;
x = 0.3*sin(t) + sin(1.3*t) + 0.9*sin(4.2*t) + 0.02*randn(1, 10001);

这是一个更不寻常的信号(均匀性不强/常规)，具有变化的频率和幅度，但通常仍为正弦波.绘制峰时，这些峰对眼睛来说显然是显而易见的，但是很难用算法来识别.

This is a more unusual (less uniform/regular) signal, with a varying frequency and amplitude, but still generally sinusoidal. The peaks are plainly obvious to the eye when plotted, but hard to identify with an algorithm.

什么是能够正确识别正弦输入信号中峰值的实时算法?在信号处理方面，我并不是真正的专家，因此获取一些信号会有所帮助.考虑正弦输入的经验法则.或者，也许我需要修改例如Jean-Paul的算法本身是为了在正弦信号上正常工作.如果是这样，将需要进行哪些修改，我将如何进行这些修改?

What is a good real-time algorithm to correctly identify the peaks in a sinusoidal input signal? I am not really an expert when it comes to signal processing, so it would be helpful to get some rules of thumb that consider sinusoidal inputs. Or, perhaps I need to modify e.g. Jean-Paul's algorithm itself in order to work properly on sinusoidal signals. If that's the case, what modifications would be required, and how would I go about making these?

推荐答案

案例1:正弦波没有噪声

如果您的正弦波不包含任何噪声，则可以使用非常经典的信号处理技术:获取一阶导数并检测其何时等于零.

Case 1: sinusoid without noise

If your sinusoid does not contain any noise, you can use a very classic signal processing technique: taking the first derivative and detecting when it is equal to zero.

例如:

function signal = derivesignal( d )

% Identify signal
signal = zeros(size(d));
for i=2:length(d)
    if d(i-1) > 0 && d(i) <= 0
        signal(i) = +1;     % peak detected
    elseif d(i-1) < 0 && d(i) >= 0
        signal(i) = -1;     % trough detected
    end
end

end

使用示例数据:

% Generate data
dt = 1/8000;
t  = (0:dt:(1-dt)/4)';
y  = sin(2*pi*60*t);

% Add some trends
y(1:1000) = y(1:1000) + 0.001*(1:1000)';
y(1001:2000) = y(1001:2000) - 0.002*(1:1000)';

% Approximate first derivative (delta y / delta x)
d = [0; diff(y)];

% Identify signal
signal = derivesignal(d);

% Plot result
figure(1); clf; set(gcf,'Position',[0 0 677 600])
subplot(4,1,1); hold on;
title('Data');
plot(t,y);
subplot(4,1,2); hold on;
title('First derivative');
area(d);
ylim([-0.05, 0.05]);
subplot(4,1,3); hold on;
title('Signal (-1 for trough, +1 for peak)');
plot(t,signal); ylim([-1.5 1.5]);
subplot(4,1,4); hold on;
title('Signals marked on data');
markers = abs(signal) > 0;
plot(t,y); scatter(t(markers),y(markers),30,'or','MarkerFaceColor','red');

这将产生:

此方法对于任何类型的正弦波都将非常有效，仅要求输入信号中不包含噪声.

This method will work extremely well for any type of sinusoid, with the only requirement that the input signal contains no noise.

一旦您的输入信号包含噪声，微分方法将失败.例如:

As soon as your input signal contains noise, the derivative method will fail. For example:

% Generate data
dt = 1/8000;
t  = (0:dt:(1-dt)/4)';
y  = sin(2*pi*60*t);

% Add some trends
y(1:1000) = y(1:1000) + 0.001*(1:1000)';
y(1001:2000) = y(1001:2000) - 0.002*(1:1000)';

% Add some noise
y = y + 0.2.*randn(2000,1);

现在将生成此结果，因为 第一差异会放大噪声. :

Will now generate this result because first differences amplify noise:

现在有很多方法可以处理噪音，最标准的方法是应用移动平均值过滤器.移动平均数的一个缺点是它们适应新信息的速度较慢，因此信号出现后可能会被识别(移动平均数有滞后性).

Now there are many ways to deal with noise, and the most standard way is to apply a moving average filter. One disadvantage of moving averages is that they are slow to adapt to new information, such that signals may be identified after they have occurred (moving averages have a lag).

另一种非常典型的方法是使用 Fourier分析确定输入数据中的所有频率，忽略所有低振幅和高频正弦波，并将其余的正弦波用作滤波器.剩余的正弦波将(大部分)从噪声中清除，然后您可以再次使用一阶微分来确定峰和谷(或者对于单个正弦波，您知道峰和谷发生在1/4和3/4 pi阶段).建议您阅读任何《信号处理理论》一书，以了解有关此技术的更多信息. Matlab对此也有一些教学材料.

Another very typical approach is to use Fourier Analysis to identify all the frequencies in your input data, disregard all low-amplitude and high-frequency sinusoids, and use the remaining sinusoid as a filter. The remaining sinusoid will be (largely) cleansed from the noise and you can then use first-differencing again to determine the peaks and troughs (or for a single sine wave you know the peaks and troughs happen at 1/4 and 3/4 pi of the phase). I suggest you pick up any signal processing theory book to learn more about this technique. Matlab also has some eductional material about this.

如果您想在硬件中使用此算法，建议您也看看WFLC(锁相环)，可以估算噪声相位wave 无进行完整的快速傅立叶变换.您可以在Wikipedia上找到用于锁相环的Matlab算法

If you want to use this algorithm in hardware, I would suggest you also take a look at WFLC (Weighted Fourier Linear Combiner) with e.g. 1 oscillator or PLL (Phase-Locked Loop) that can estimate the phase of a noisy wave without doing a full Fast Fourier Transform. You can find a Matlab algorithm for a phase-locked loop on Wikipedia.

在这里，我建议采用一种稍微复杂些的方法，该方法可以实时识别峰值和谷值:通过移动

I will suggest a slightly more sophisticated approach here that will identify the peaks and troughs in real-time: fitting a sine wave function to your data using moving least squares minimization with initial estimates from Fourier analysis.

这是我的职责:

function [result, peaks, troughs] = fitsine(y, t, eps)

% Fast fourier-transform
f = fft(y);
l = length(y);
p2 = abs(f/l);
p1 = p2(1:ceil(l/2+1));
p1(2:end-1) = 2*p1(2:end-1);
freq = (1/mean(diff(t)))*(0:ceil(l/2))/l;

% Find maximum amplitude and frequency
maxPeak = p1 == max(p1(2:end)); % disregard 0 frequency!
maxAmplitude = p1(maxPeak);     % find maximum amplitude
maxFrequency = freq(maxPeak);   % find maximum frequency

% Initialize guesses
p = [];
p(1) = mean(y);         % vertical shift
p(2) = maxAmplitude;    % amplitude estimate
p(3) = maxFrequency;    % phase estimate
p(4) = 0;               % phase shift (no guess)
p(5) = 0;               % trend (no guess)

% Create model
f = @(p) p(1) + p(2)*sin( p(3)*2*pi*t+p(4) ) + p(5)*t;
ferror = @(p) sum((f(p) - y).^2);
% Nonlinear least squares
% If you have the Optimization toolbox, use [lsqcurvefit] instead!
options = optimset('MaxFunEvals',50000,'MaxIter',50000,'TolFun',1e-25);
[param,fval,exitflag,output] = fminsearch(ferror,p,options);

% Calculate result
result = f(param);

% Find peaks
peaks = abs(sin(param(3)*2*pi*t+param(4)) - 1) < eps;

% Find troughs
troughs = abs(sin(param(3)*2*pi*t+param(4)) + 1) < eps;

end

如您所见，我首先执行傅里叶变换来找到数据的幅度和频率.然后，我使用模型 a + b sin(ct + d)+ et 将正弦曲线拟合到数据中.拟合值表示一个正弦波，我知道+1和-1分别是峰值和谷值.因此，我可以将这些值标识为信号.

As you can see, I first perform a Fourier transform to find initial estimates of the amplitude and frequency of the data. I then fit a sinusoid to the data using the model a + b sin(ct + d) + et. The fitted values represent a sine wave of which I know that +1 and -1 are the peaks and troughs, respectively. I can therefore identify these values as the signals.

这对于具有(缓慢变化的)趋势和一般(白色)噪声的正弦曲线非常有用:

This works very well for sinusoids with (slowly changing) trends and general (white) noise:

% Generate data
dt = 1/8000;
t  = (0:dt:(1-dt)/4)';
y  = sin(2*pi*60*t);

% Add some trends
y(1:1000) = y(1:1000) + 0.001*(1:1000)';
y(1001:2000) = y(1001:2000) - 0.002*(1:1000)';

% Add some noise
y = y + 0.2.*randn(2000,1);

% Loop through data (moving window) and fit sine wave
window = 250;   % How many data points to consider
interval = 10;  % How often to estimate
result = nan(size(y));
signal = zeros(size(y));
for i = window+1:interval:length(y)
    data = y(i-window:i);   % Get data window
    period = t(i-window:i); % Get time window
    [output, peaks, troughs] = fitsine(data,period,0.01);

    result(i-interval:i) = output(end-interval:end);
    signal(i-interval:i) = peaks(end-interval:end) - troughs(end-interval:end);
end

% Plot result
figure(1); clf; set(gcf,'Position',[0 0 677 600])
subplot(4,1,1); hold on;
title('Data');
plot(t,y); xlim([0 max(t)]); ylim([-4 4]);
subplot(4,1,2); hold on;
title('Model fit');
plot(t,result,'-k'); xlim([0 max(t)]); ylim([-4 4]);
subplot(4,1,3); hold on;
title('Signal (-1 for trough, +1 for peak)');
plot(t,signal,'r','LineWidth',2); ylim([-1.5 1.5]);
subplot(4,1,4); hold on;
title('Signals marked on data');
markers = abs(signal) > 0;
plot(t,y,'-','Color',[0.1 0.1 0.1]);
scatter(t(markers),result(markers),30,'or','MarkerFaceColor','red');
xlim([0 max(t)]); ylim([-4 4]);

这种方法的主要优点是:

Main advantages of this approach are:

• 您具有数据的实际模型，因此您可以在信号发生之前预测未来！ (例如，修复模型并通过输入将来的时间段来计算结果)
• 您不需要每个周期都估算模型(请参见代码中的参数interval)

缺点是您需要选择回溯window，但是任何用于实时检测的方法都会遇到此问题.

The disadvantage is that you need to select a lookback window, but you will have this problem with any method that you use for real-time detection.

Data是输入数据，Model fit是数据的拟合正弦波(请参见代码)，Signal指示波峰和波谷，而Signals marked on data则给人一种算法精确度的印象. 注意:观看模型拟合以使其自身适应图中间的趋势！

Data is the input data, Model fit is the fitted sine wave to the data (see code), Signal indicates the peaks and troughs and Signals marked on data gives an impression of how accurate the algorithm is. Note: watch the model fit adjust itself to the trend in the middle of the graph!

那应该让您入门.也有很多关于信号检测理论的优秀书籍(仅用Google来形容)，将深入探讨这些类型的技术.祝你好运！

That should get you started. There are also a lot of excellent books on signal detection theory (just google that term), which will go much further into these types of techniques. Good luck!

这篇关于嘈杂正弦时间序列中的实时峰值检测的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！