function [ampl,phase]=wave_fit(t,y,period)
% Аппроксимаци? р?да набором синусоид с заданными периодами
%
%            [ampl,phase]=wave_fit(t,y,period)
%
%     t - вектор абсцисс
%     y - вектор измеренных значений
%period - вектор периодов аппроксимирующих синусоид.
%         Его длина определ?ет кол-во гармоник, единица
%         измерени? - та же, что и дл? t.
%
%  ampl - вектор амплитуд аппроксимирующих гармоник
% phase - вектор фаз (в тех же единицах, что и t)
%
%!!! дл? v3.5 произведение length(t) * (2*length(period)) д.б. < 8188 !!!
%
%См. также:   WAVE_VAL

%
% Система SUM{ c(k)*SIN[f(k)*t+ph(k)] } линеаризируетс? к виду
% SUM{ a(k)*SIN[f(k)*t] + b(k)*COS[f(k)*t]}, где f(k) - частота
% k-ой гармоники, a(k)=c(k)*COS(ph(k)), b(k)=c(k)*SIN(ph(k)).
% Решение A*X=Y ищетс? с помощью МНК: X=A\Y.

%
% Хомутов С.Ю.,   12.04.96
%                 05.11.97 Ревизи? дл? v.4.2
%
freq=2*pi./period;      % частоты гармоник
nwave=length(period);   % кол-во гармоник дл? аппроксимации
if sscanf(version,'%d',1)==3 & length(t)*2*nwave>8188,
   error('Не смогу создать матрицу коэффициентов дл? ver.3')
end
% Инициализаци? матрицы коэффициентов
a=zeros(length(t),2*nwave);
for k=1:nwave
   a(:,2*k-1)=sin(freq(k)*t);
   a(:,2*k)  =cos(freq(k)*t);
end
% Решение системы a*x=y ---> x=a\y методом наименьших квадратов
x=a\y;
for k=1:nwave
   ampl(k)=sqrt(x(2*k-1)^2+x(2*k)^2);
   phase(k)=atan2(x(2*k),x(2*k-1))/freq(k);
end
