function [X,Xsd,P,Psd,xac,xdc,res,tw] = osc(x,t,f,w,o,tw)

% OSC   Modulation analysis
%
%   [XF,XFSD,P,PSD,XAC,XDC,RES,TW] = OSC(X,T,F[,W,O,TW]) fits to the time signal
%   X = XF*exp(-2*pi*j*F*T) + P*XP = XAC + XDC
%
%   X    signal
%   T    signal time base
%   F    arrays of frequencies
%   W    width of analysis window, default or [] one period of the smallest F
%   O    polynomial term orders, default or [] no polynomial
%        or [-O TO] for B-splines with order O and knot separation TO
%   TW   desired times, if absent takes T, [] gives a window spacing
%
%   XF   complex amplitude at each (TW,F,column of X)
%   XFSD standard deviation of complex amplitudes
%   P    polynomial coefficients at each (TW,O,column of X)
%   PSD  standard deviation of polynomial coefficients
%   XAC  fitted ac signal
%   XDC  fitted dc signal
%   RES  residual after fitting
%   TW   fitting time
%
% [+GenLib General Purpose Library+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% demo
if nargin == 0
 t = linspace(0,1,101)';
 f = [5 6];
 x = [cos(2*pi*t*f) sin(2*pi*t*f)] * [1 0;0 0;0 0;0 1] + t*[1 2] + 1 + randn(101,2)/10;
 [X,Xsd,P,Psd,xac,xdc,res,tw] = osc(x,t,f,.5,[0 1],linspace(.25,.75,21));
 plot(t,x,tw,xac+xdc,tw,res)
 return
end

% missing argumentws
if nargin < 4,  w = 1/min(f); end
if isempty(x),  w = 1/min(f); end
if nargin < 5,  o = []; end
if nargin < 6,  tw = t; end
if isempty(tw), tw = t(1)+w/2:w:t(length(t))-w/2; end

% align array if necessary
t  = t(:);
tw = tw(:);
f  = f(:)';

% spline base line
if any(o < 0)
 nb = -o(1);
 tb = (-nb:floor((t(end)-(t(1)))/o(2))+nb)*o(2)+t(1);
 no = length(tb)-nb;
else
 no = length(o);
end

% fuq
[nt,nx] = size(x);
nf  = length(f);
nv  = 2*nf + no;
ntw = length(tw);
w   = w/2;

% dependant variable matrix
v = zeros(nt,nv);
v(:,1:nf)   = t*(2*pi*f);
v(:,1:2*nf) = [cos(v(:,1:nf)) sin(v(:,1:nf))];
if any(o < 0)
 v(:,2*nf+1:end) = bspbase(tb,nb,t)';
else
 for k = 1:no, v(:,2*nf+k) = t.^o(k); end
end

% create dependant variables matrix at wanted times tw
vw = zeros(ntw,nv);
vw(:,1:nf)   = tw*(2*pi*f);
vw(:,1:2*nf) = [cos(vw(:,1:nf)) sin(vw(:,1:nf))];
if any(o < 0)
 vw(:,2*nf+1:end) = bspbase(tb,nb,tw)';
else
 for k = 1:no, v(:,2*nf+k) = t.^o(k); end
end

% preallocation
X   = repmat(NaN,[ntw,nf,nx]);
Xsd = repmat(NaN,[ntw,nf,nx]);
P   = repmat(NaN,[ntw,no,nx]);
Psd = repmat(NaN,[ntw,no,nx]);
xac = repmat(NaN,[ntw,nx]);
xdc = repmat(NaN,[ntw,nx]);
res = repmat(NaN,[ntw,nx]);

% indices for distribution of the results
kr = 1:nf;
ki = kr+nf;
ko = 2*nf+(1:no);

% scan the analysis window
for k = 1:ntw

 % window sample selection
 sub   = find(abs(t-tw(k)) <= w);
 subv = v(sub,:);
 kv = find(any(subv ~= 0));
 
 if length(sub) >= length(kv) % enough samples for parameter number
  subx = x(sub,:);

  % computes coefficients and invert covariance matrix
  b = zeros(nv,nx); c = zeros(nv,nv);
  b(kv,:)  = xmldivide(subv(:,kv),subx);
  c(kv,kv) = inv(subv(:,kv)'*subv(:,kv));

  % fitted value and residual after fitting
  xac(k,:) = vw(k,1:2*nf)    * b(1:2*nf,:);
  xdc(k,:) = vw(k,2*nf+1:nv) * b(2*nf+1:nv,:);
  res(k,:) = sqrt(sum((subv*b-subx).^2)/(length(sub)-nv));
  bsd = sqrt(diag(c)) * res(k,:);

  % sort out the results
  X  (k,:,:) = b  (kr,:) + j*b  (ki,:);
  Xsd(k,:,:) = bsd(kr,:) + j*bsd(ki,:);
  P  (k,:,:) = b  (ko,:);
  Psd(k,:,:) = bsd(ko,:);
 end

end

