Matlab：帮助理解正弦曲线拟合

Question

我有一个未知的正弦波，其中包含一些我正在尝试重建的噪声。最终目标是提出一个 C 算法来找到正弦波的幅度、直流偏移、相位和频率，但我首先在 Matlab（实际上是 Octave）中进行原型设计。正弦波的形式是

y = a + b*sin(c + 2*pi*d*t)
a = dc offset
b = amplitude
c = phase shift (rad)
d = frequency

我找到了this example，并且在 cmets 中，John D'Errico 提出了一种使用最小二乘法将正弦波拟合到数据的方法。这是一个简洁的小算法，效果非常好，但我很难理解一个方面。算法如下：

算法

假设您有如下形式的正弦波：

(1) y = a + b*sin(c+d*x)

使用身份

(2) sin(u+v) = sin(u)*cos(v) + cos(u)*sin(v)

我们可以将(1)改写为

(3) y = a + b*sin(c)*cos(d*x) + b*cos(c)*sin(d*x)

由于b*sin(c) 和b*cos(c) 是常量，因此可以将它们包装成常量b1 和b2。

(4) y = a + b1*cos(d*x) + b2*sin(d*x)

这是用于拟合正弦波的方程。创建一个函数来生成回归系数和平方和残差。

(5) cfun = @(d) [ones(size(x)), sin(d*x), cos(d*x)] \ y;
(6) sumerr2 = @(d) sum((y - [ones(size(x)), sin(d*x), cos(d*x)] * cfun(d)) .^ 2);

接下来，sumerr2 使用 fminbnd 最小化频率 d，下限 l1 和上限 l2。

(7) dopt = fminbnd(sumerr2, l1, l2);

现在可以计算 a、b 和 c。计算 a、b 和 c 的系数由 (4) 在 dopt 处给出

(8) abb = cfun(dopt);

直流偏移只是第一个值

(9) a = abb(1);

使用三角标识来查找b

(10) sin(u)^2 + cos(u)^2 = 1
(11) b = sqrt(b1^2 + b2^2)
(12) b = norm(abb([2 3]));

终于找到了相位偏移

(13) b1 = b*cos(c)
(14) c = acos(b1 / b);
(15) c = acos(abb(2) / b);

问题

(5) 和 (6) 中发生了什么？有人可以分解伪代码中发生的事情，或者以更明确的方式执行相同的功能吗？

(5) cfun = @(d) [ones(size(x)), sin(d*x), cos(d*x)] \ y;
(6) sumerr2 = @(d) sum((y - [ones(size(x)), sin(d*x), cos(d*x)] * cfun(d)) .^ 2);

另外，鉴于（4）不应该是：

[ones(size(x)), cos(d*x), sin(d*x)]

代码

这是完整的 Matlab 代码。蓝线是实际信号。绿线是重构信号。

close all
clear all

y  = [111,140,172,207,243,283,319,350,383,414,443,463,483,497,505,508,503,495,479,463,439,412,381,347,311,275,241,206,168,136,108,83,63,54,45,43,41,45,51,63,87,109,137,168,204,239,279,317,348,382,412,439,463,479,496,505,508,505,495,483,463,441,414,383,350,314,278,245,209,175,140,140,110,85,63,51,45,41,41,44,49,63,82,105,135,166,200,236,277,313,345,379,409,438,463,479,495,503,508,503,498,485,467,444,415,383,351,318,281,247,211,174,141,111,87,67,52,45,42,41,45,50,62,79,104,131,163,199,233,273,310,345,377,407,435,460,479,494,503,508,505,499,486,467,445,419,387,355,319,284,249,215,177,143,113,87,67,55,46,43,41,44,48,63,79,102,127,159,191,232,271,307,343,373,404,437,457,478,492,503,508,505,499,488,470,447,420,391,360,323,287,254,215,182,147,116,92,70,55,46,43,42,43,49,60,76,99,127,159,191,227,268,303,339,371,401,431,456,476,492,502,507,507,500,488,471,447,424,392,361,326,287,287,255,220,185,149,119,92,72,55,47,42,41,43,47,57,76,95,124,156,189,223,258,302,337,367,399,428,456,476,492,502,508,508,501,489,471,451,425,396,364,328,294,259,223,188,151,119,95,72,57,46,43,44,43,47,57,73,95,124,153,187,222,255,297,335,366,398,426,451,471,494,502,507,508,502,489,474,453,428,398,367,332,296,262,227,191,154,124,95,75,60,47,43,41,41,46,55,72,94,119,150,183,215,255,295,331,361,396,424,447,471,489,500,508,508,502,492,475,454,430,401,369,335,299,265,228,191,157,126,99,76,59,49,44,41,41,46,55,72,92,118,147,179,215,252,291,328,360,392,422,447,471,488,499,507,508,503,493,477,456,431,403]';
fs = 100e3;
N  = length(y);
t  = (0:1/fs:N/fs-1/fs)';

cfun    = @(d) [ones(size(t)), sin(2*pi*d*t), cos(2*pi*d*t)]\y;
sumerr2 = @(d) sum((y - [ones(size(t)), sin(2*pi*d*t), cos(2*pi*d*t)] * cfun(d)) .^ 2);
dopt    = fminbnd(sumerr2, 2300, 2500);
abb     = cfun(dopt);

a = abb(1);
b = norm(abb([2 3]));
c = acos(abb(2) / b);
d = dopt;

y_reconstructed = a + b*sin(2*pi*d*t - c);

figure(1)
hold on
title('Signal Reconstruction')
grid on
plot(t*1000, y, 'b')
plot(t*1000, y_reconstructed, 'g')
ylim = get(gca, 'ylim');
xlim = get(gca, 'xlim');
text(xlim(1), ylim(2) - 15, [num2str(b) ' cos(2\pi * ' num2str(d) 't - ' ... 
                             num2str(c * 180/pi) ') + ' num2str(a)]);
hold off

Answer 1

(5) 和 (6) 定义了可以在优化代码中使用的匿名函数。 cfun 返回一个数组，它是 t、y 和参数 d（即会变化的优化参数）的函数。同样， sumerr2 是另一个匿名函数，具有相同的参数，这次返回一个标量。该标量将是 fminbnd 最小化的误差。