16

我看过这个和这个。

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

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

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

编辑：我的曲线真的很吵，所以用 FFT 来计算频率可能会导致误报和不合适的结果。我正在寻找一种更稳健的拟合方式。

约 170 点的曲线

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

在此处输入图像描述

噪声明显是非高斯的，并且与数据的幅度相比很大。我已经在高斯分布的噪声上尝试过 FFT，我的配合非常好。在这里，它的失败非常严重。

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

4

4 回答 4

7

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

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

更新：

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

第 11 列中的数据图 第 13 列中的数据图

我使用 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 是否可以从噪声信号中正确估计值。该信号由两条添加的正弦曲线组成，频率分别为 50 Hz、120 Hz，幅度分别为 0.7、1.0。下图中，红色曲线是原始曲线，蓝色曲线是添加了噪点的曲线。（图被裁剪）

PHD精度测试

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 估计 amp/freq 值为：

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)

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

回复@亚历克斯：

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

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

于 2013-01-24T17:11:25.670 回答
4

Kitchi：你能提供一些样本数据吗？您必须使用的典型信号多长时间？（就样本数和正弦波周期数而言） 信噪比以 dB 为单位是多少？

  1. 在你知道什么算法会起作用之前，我建议你在 python/numpy/scipy（或 matlab/octave，或 R/S，或 Mathematica ......）中制作原型，无论你最喜欢什么原型语言/工具集，除了 C。它将节省大量时间，并且您将使用更丰富的工具。

  2. 您确定噪声会严重影响 FFT 吗？这不一定是一个好的假设，特别是如果噪声相对“白”，并且分析窗口很长。如果正弦波的频率非常稳定，你可以做一个巨大的 FFT 并将信号从比信号强几个数量级的噪声中提取出来。尝试预期正弦波的几百到几百万个周期。

  3. 曲线拟合正弦波效果不佳。我猜周期性会产生很多局部最小值，而相移变量也会使问题显着非线性。您可以从下面链接的遇到相同问题的其他人那里看到一些问题。你最好尝试几乎任何其他东西而不是非线性最小二乘拟合，除非你预先线性化问题，这让我......

  4. 自相关非常适合这种事情。尝试一次计算整个信号的自相关（如果源频率稳定，数据越多越好）。正弦波周期作为自相关中的一个高峰应该非常明显，并且您可能会比使用 FFT 获得更准确的频率估计（除非您使用非常大的 FFT）。此外，您可以从第一个高自相关峰的高度计算平均幅度。

编辑：经过进一步研究，有更多的技术可能比 FFT 更适合您的问题。Pisarenko 的调和分解（下面 Fredrik Rubin 首次提出）是一个；另一个是最小二乘谱分析（LSSA），它与您最初的问题想法非常相似。LSSA 有许多变体，例如 Lomb-Scargle、基础追踪等，它们以各种方式处理我上面描述的拟合问题。但是我认为，如果您在大型 FFT 中绝对看不到任何信号，那么其他任何方法都不太可能找到任何东西:)

PS 有关无法很好地拟合正弦波的其他问题，请参阅：

于 2013-01-21T05:23:10.967 回答
2

如果您要对 sin 进行回归，则可以使用 FFT 应用傅立叶变换。

编辑

尝试使用过滤器去除噪音。如果您有传感器等物理源，请在传感器上放置低通滤波器。FFT 是相对较差的滤波器。

EDIT2 - 这个测量是完全错误的

可能是您进行了错误的测量。根据Nyquist-Shannon 采样定理，您的采样频率太低，或者输入频率太高。这会导致错误的解决方案，因为如果您使用 5kHz 采样进行采样，例如 3kHz，您将根据该定理测量 2kHz。

我确信您无法通过这种测量来判断正确的输入频率。

于 2013-01-07T14:21:56.363 回答
2

这实际上是一个谱估计问题。您正在尝试估计一个“线谱”，您知道您拥有的正弦波数量（在您的情况下是一个）。像MUSIC或ESPRIT这样的方法应该能够解决这个问题。

作为参考，Stoica的书会派上用场。本书的第 4 章是线谱的参数方法，其中包含用于查找所需信号的幅度、相位和频率的算法。这本书还附带了在 MATLAB 中实现的算法，它们也很容易自己实现。

于 2013-01-27T13:59:58.760 回答