c - Matlab: Help understanding sinusoidal curve fit -

i have unknown sine wave noise trying reconstruct. ultimate goal come c algorithm find amplitude, dc offset, phase, , frequency of sine wave prototyping in matlab (octave actually) first. sine wave of form

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

i have found this example , in comments john d'errico presents method using least squares fit sine wave data. neat little algorithm , works remarkably having difficulties understanding 1 aspect. algorithm follows:

algorithm

suppose have sine wave of form:

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

using identity

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

we can rewrite (1) as

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

since b*sin(c) , b*cos(c) constants, these can wrapped constants b1 , b2.

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

this equation used fit sine wave. function created generate regression coefficients , 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 minimized frequency d using fminbnd lower limit l1 , upper limit l2.

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

now a, b, , c can computed. coefficients compute a, b, , c given (4) @ dopt

(8) abb = cfun(dopt); 

the dc offset first value

(9) = abb(1); 

a trig identity used find b

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

finally phase offset found

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

question

what going on in (5) , (6)? can break down happening in pseudo-code or perhaps perform same function in 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 be:

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

code

here matlab code in full. blue line actual signal. green line reconstructed signal.

close clear  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);  = abb(1); b = norm(abb([2 3])); c = acos(abb(2) / b); d = dopt;  y_reconstructed = + 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 

(5) , (6) defining anonymous functions can used within optimisation code. cfun returns array function of t, y , parameter d (that optimisation parameter varied). similarly, sumerr2 anonymous function, same arguments, time returning scalar. scalar error minimised fminbnd.