c - 正弦曲线的线性/非线性拟合

Question

我看过this和 this .

但我有一个稍微不同的问题。我知道我的数据是一条正弦曲线，周期未知，振幅未知，具有附加的非高斯分布噪声。

我正在尝试使用 GSL 来适应它C 中的非线性算法，但拟合绝对糟糕。我想知道我是否(错误地)使用了非线性拟合算法，而我应该使用线性拟合算法？

如何判断特定数据集需要线性算法还是非线性算法？

编辑:我的曲线确实有噪声，因此采用 FFT 计算频率可能会导致误报和不合适。我正在寻找一种稍微更稳健的拟合方式。

如您所见，上面的图大约有 170 个点，下面的图大约有 790 个点。

噪声明显是非高斯分布的，并且与数据的振幅相比很大。我试过对高斯分布噪声进行 FFT，我的拟合度非常好。在这里，它失败得很严重。

添加链接到first时间序列数据。文件中的每一列都是不同的时间序列。

Answer 1

如果您知道您的数据是一条正弦曲线(可以表示为多个复指数)，那么您可以使用 Pisarkenko 的调和分解； http://en.wikipedia.org/wiki/Pisarenko_harmonic_decomposition

但是，如果您可以访问更多数据点，我的方法仍然会使用 DFT。

更新:

我对您的数据使用了 Pisarenko 的谐波分解 (PHD)，即使您的信号非常短(每个只有 86 个数据点)，但如果有更多可用数据，PHD 算法肯定具有潜力。下面包括 24 个信号中的两个(数据的第 11 和 13 列)，用蓝色表示，红色的正弦曲线对应于 PHD 的估计幅度/频率值。 (注意相移是未知的)

我使用 MATLAB (pisar.m) 执行 PHD:http://www.mathworks.com/matlabcentral/fileexchange/74

% assume data is one single sine curve (in noise)
SIN_NUM = 1; 

for DATA_COLUMN = 1:24
    % obtain amplitude (A), and frequency (f = w/2*pi) estimate
    [A f]=pisar(data(:,DATA_COLUMN),SIN_NUM);

    % recreated signal from A, f estimate
    t = 0:length(data(:,DATA_COLUMN))-1;
    y = A*cos(2*pi*f*t);

    % plot original/recreated signal
    figure; plot(data(:,DATA_COLUMN)); hold on; plot(y,'r')
    title({'data column ',num2str(DATA_COLUMN)});

    disp(A)
    disp(f)
end

这导致了

1.9727     % amp. for  column 11
0.1323     % freq. for column 11
2.3231     % amp. for  column 13
0.1641     % freq. for column 13

PHD 验证:

我还做了另一个测试，我知道振幅和频率的值，然后添加噪声以查看 PHD 是否可以从噪声信号中正确估计值。该信号由两条附加的正弦曲线组成，频率分别为 50 Hz、120 Hz，振幅分别为 0.7、1.0。下图中，红色曲线为原始曲线，蓝色曲线为添加噪声后的曲线。 (图被裁剪)

Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector

% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 0.4*randn(size(t)); % Sinusoids plus noise

figure;
plot(Fs*t(1:100),y(1:100)); hold on; plot(Fs*t(1:100),x(1:100),'r')
title('Signal Corrupted with Zero-Mean Random Noise (Blue), Original (Red)')

[A, f] = pisar(y',2); 
disp(A)
disp(f/Fs)

PHD 估计放大器/频率值为:

0.7493    % amp wave 1  (actual 0.7)
0.9257    % amp wave 2  (actual 1.0)
58.5      % freq wave 1 (actual 50)
123.8     % freq wave 2 (actual 120)

对于相当多的噪音来说还不错，而且只知道信号包含的波数。

回复@Alex:

是的，这是一个很好的算法，我在 DSP 研究期间遇到过它，我认为它工作得很好，但重要的是要注意 Pisarenko 的 Harm.Dec。将任何信号建模为 N > 0 正弦波，N 从一开始就指定，并使用该值忽略噪声。因此，根据定义，只有当您粗略地知道您的数据由正弦曲线组成时，它才有用。如果您对 N 的值一无所知，并且需要针对一千个不同的值运行该算法，那么绝对推荐使用不同的方法。也就是说，此后的评估很简单，因为它返回 N 个振幅和频率值。

多信号分类 (MUSIC) 是另一种算法，它在 Pisarenko 中断的地方继续进行。 http://en.wikipedia.org/wiki/Multiple_signal_classification