Matlab: Help understanding sinusoidal curve fit

Question

I have an unknown sine wave with some noise that I am trying to reconstruct. The ultimate goal is to come up with a C algorithm to find the amplitude, dc offset, phase, and frequency of a sine wave but I am prototyping in Matlab (Octave actually) first. The sine wave is of the form

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

I have found this example and in the comments John D'Errico presents a method for using Least Squares to fit a sine wave to data. It is a neat little algorithm and works remarkably well but I am having difficulties understanding one aspect. The algorithm is as follows:

Algorithm

Suppose you have a sine wave of the form:

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

Using the identity

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

We can rewrite (1) as

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

Since b*sin(c) and b*cos(c) are constants, these can be wrapped into constants b1 and b2.

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

This is the equation that is used to fit the sine wave. A function is created to generate regression coefficients and a sum-of-squares residual error.

(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);

Next, sumerr2 is minimized for the frequency d using fminbnd with lower limit l1 and upper limit l2.

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

Now a, b, and c can be computed. The coefficients to compute a, b, and c are given from (4) at dopt

(8) abb = cfun(dopt);

The dc offset is simply the first value

(9) a = abb(1);

A trig identity is used to find b

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

Finally the phase offset is found

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

Question

What is going on in (5) and (6)? Can someone break down what is happening in pseudo-code or perhaps perform the same function in a more explicit way?

(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);

Also, given (4) shouldn't it be:

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

Code

Here is the Matlab code in full. Blue line is the actual signal. Green line is the reconstructed signal.

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) and (6) are defining anonymous functions that can be used within the optimisation code. cfun returns an array that is a function of t, y and the parameter d (that is the optimisation parameter that will be varied). Similarly, sumerr2 is another anonymous function, with the same arguments, this time returning a scalar. That scalar will be the error that is to be minimised by fminbnd.