案例1：正弦波没有噪声
如果您的正弦曲线不包含任何噪声，您可以使用非常经典的信号处理技术：取一阶导数并检测它何时等于零。
例如：
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');
这会产生：
这种方法对于任何类型的正弦曲线都非常有效，唯一的要求是输入信号不包含噪声。
案例2：正弦有噪声
一旦您的输入信号包含噪声，导数方法就会失败。例如：
% 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);
现在会生成这个结果，因为first differences amplify noise:
现在处理噪音的方法有很多，最标准的方法是申请moving average filter。移动平均线的一个缺点是它们适应新信息的速度很慢，因此信号可能会在它们发生后被识别（移动平均线有滞后）。
另一种非常典型的方法是使用Fourier Analysis 来识别输入数据中的所有频率，忽略所有低幅度和高频正弦波，并使用剩余的正弦波作为滤波器。剩余的正弦波将（大部分）从噪声中清除，然后您可以再次使用一阶差分来确定波峰和波谷（或者对于单个正弦波，您知道波峰和波谷发生在 1/4 和 3/4 pi相）。我建议您阅读任何信号处理理论书籍以了解有关此技术的更多信息。 Matlab 也有一些关于这个的educational material。
如果你想在硬件中使用这个算法，我建议你也看看 WFLC (Weighted Fourier Linear Combiner)，例如1 个振荡器或 PLL (Phase-Locked Loop)， 无需进行完整的快速傅立叶变换即可估计噪声波的相位。您可以找到锁相环的 Matlab 算法on Wikipedia。
我将在这里建议一种稍微复杂一点的方法，该方法将实时识别波峰和波谷：使用移动 least squares minimization 和傅立叶分析的初步估计，为您的数据拟合正弦波函数。
这是我的功能：
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
如您所见，我首先执行Fourier transform 来查找数据幅度和频率的初始估计值。然后，我使用模型 a + b sin(ct + d) + et 将正弦曲线拟合到数据中。拟合值代表一个正弦波，我知道 +1 和 -1 分别是波峰和波谷。因此，我可以将这些值识别为信号。
这对于具有（缓慢变化）趋势和一般（白）噪声的正弦曲线非常有效：
% 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]);
这种方法的主要优点是：
- 您拥有数据的实际模型，因此您可以在信号发生之前预测未来的信号！ （例如修复模型并通过输入未来时间段来计算结果）
- 您不需要每期都估计模型（见代码中的参数
interval）
缺点是需要选择回溯window，但是任何用于实时检测的方法都会有这个问题。
视频演示
Data 是输入数据，Model fit 是拟合到数据的正弦波（见代码），Signal 表示波峰和波谷，Signals marked on data 表示算法的准确度。 注意：观察模型拟合调整到图表中间的趋势！
这应该可以帮助您入门。还有很多关于信号检测理论的优秀书籍（只需 google 那个术语），它们将更深入地介绍这些类型的技术。祝你好运！