拟合具有失真时基的正弦波

Question

我想知道在 Matlab 中拟合具有失真时基的正弦波的最佳方法。

时间失真由 n 阶多项式 (n~10) 给出，形式为t_distort = P(t)。

例如，考虑失真t_distort = 8 + 12t + 6t^2 + t^3（这只是(t-2)^3的幂级数展开式）。

这将使正弦波失真如下：

我希望能够找到这个失真正弦波的失真。 （即我想找到函数t = G(t_distort)，但是t_distort = P(t)是未知的。）

Answer 1

如果你的分辨率足够高，那么这基本上是一个角度解调问题。解调角度调制信号的标准方法是求导数，然后是包络检波器，然后是积分器。

由于我不知道您使用的确切数字，我将用我自己的数字举一个例子。假设我原来的时基从 0 到 100 有 1000 万个点：

t = 0:0.00001:100;

然后我得到失真的时基并计算失真的正弦波：

td = 0.02*(t+2).^3;
yd = sin(td);

现在我可以解调它了。使用近似差除以之前的步长来取“导数”：

ydot = diff(yd)/0.00001;

信封可以是easily detected:

envelope = abs(hilbert(ydot));

这给出了 P(t) 导数的近似值。最后一步是积分器，我可以使用累积和来近似（我们必须再次按步长对其进行缩放）：

tdguess = cumsum(envelope)*0.00001;

这给出了一条与原始失真时基几乎相同的曲线（因此，它给出了 P(t) 的良好近似值）：

您将无法获得多项式的常数项，因为我们从它的导数中进行了近似，这当然消除了常数项。您甚至无法仅从 yd 在数学上找到唯一的常数项，因为无限多的值将产生相同的失真正弦波。如果您知道 P(t) 的度数，则可以使用 polyfit 获得 P(t) 的其他三个系数（忽略最后一个数字，它是常数项）：

>> polyfit(t(1:10000000), tdguess, 3)

ans =

    0.0200    0.1201    0.2358    0.4915

这非常接近原始值，除了常数项：0.02*(t+2)^3 = 0.02t^3 + 0.12t^2 + 0.24t + 0.16。

你想要逆映射 Q(t)。你能知道目前发现的 P(t) 的近似值吗？

Answer 2

这是一个分析驱动的路径，它采用信号的asin 并正确展开角度。然后您可以在角度上使用polyfit 或使用其他拟合方法拟合多项式（搜索fit 并查看）。最后，对拟合函数取一个罪，并将信号与拟合函数进行比较......请参阅这个教学示例：

% generate data
t=linspace(0,10,1e2);
x=0.02*(t+2).^3;
y=sin(x);

% take asin^2 to obtain points of "discontinuity" where then asin hits +-1
da=(asin(y).^2);
[val locs]=findpeaks(da); % this can be done in many other ways too...

% construct the asin according to the proper phase unwrapping
an=NaN(size(y));
an(1:locs(1)-1)=asin(y(1:locs(1)-1));
for n=2:numel(locs)
    an(locs(n-1)+1:locs(n)-1)=(n-1)*pi+(-1)^(n-1)*asin(y(locs(n-1)+1:locs(n)-1));
end
an(locs(n)+1:end)=n*pi+(-1)^(n)*asin(y(locs(n)+1:end));

r=~isnan(an);
p=polyfit(t(r),an(r),3);

figure;  
subplot(2,1,1); plot(t,y,'.',t,sin(polyval(p,t)),'r-');
subplot(2,1,2); plot(t,x,'.',t,(polyval(p,t)),'r-');
title(['mean error ' num2str(mean(abs(x-polyval(p,t))))]);

p =

    0.0200    0.1200    0.2400    0.1600

我使用NaN 预先分配并避免在不连续点（位置）使用asin 的原因是为了减少以后拟合的错误。如您所见，对于 0,10 之间的 100 个点，平均误差是浮点精度的数量级，并且多项式系数尽可能精确。

您不采用导数这一事实（如在非常优雅的希尔伯特变换中）允许在数值上精确。在相同条件下，希尔伯特变换解决方案的平均误差会大得多（单位阶与 1e-15 的阶数）。

此方法的唯一限制是您需要在 asin 翻转方向的机制中有足够的点，并且 sin 内的函数表现良好。如果存在采样问题，您可以截断数据并仅保持更接近零的范围，这样就足以表征sin 中的函数。毕竟，您不需要数百万个操作点来适应 3 参数函数。