Разложить спектр на индивидуальные гауссовы полосы

function [M,amplitudes]=cal_ampl(y,m,std,figure_number)
%расчет амплитуд по методу наименьших квадратов
%y - массив измеренных данных
%m - массив местоположений(в местоположениях максимума стоит 1)
%std - массив стандартных отклонений(в местоположении стоит стандартное отклонение)
%figure_number -номер фигуры, если много раз делать, мне кажется так удобнее
n=length(y);
M=0;%сначала количество гауссиан нулевое
amplitudes=zeros(1,1);%сначала массив амплитуд нулевой
for i=1:n,
    if (m(i)==1), M=M+1;end;%в этой точке есть гауссиана
end;
%-------------------------------------------------------------
if (M>0),%если есть хотя бы одна гауссиана - иначе выход с нулевыми значениями
    mo=zeros(M,1);%рабочий массив местоположений гауссиан
    std=zeros(M,M);%рабочий массив полуширин гауссиан
    %заносим значения в рабочие массивы
    m=1;
    for i=1:n,
      if (ar_m(i)==1), 
      mo(m)=i;sigma(m)=ar_sigma(i);    
      m=m+1;
      end;
    end;
    mat=zeros(M,M);%матрица для обращения
    B=zeros(M,1);%столбец свободных членов
    phi=zeros(M,n);%массив гауссиан
    %формируем массив гауссиан
    for j=1:M
        for i=1:n,
            phi(j,i)=exp(-(i-mo(j))^2/(2*sigma(j)^2));%формула гаусстаны
        end;
    end;
    %формируем матрицу для обращения и столбец свободных членов
    for j=1:M
        s=0;for i=1:n,s=s+g(i)*phi(j,i);end;B(j)=s;
        for k=1:M
            mat(j,k)=sum(phi(j,:).*phi(k,:));
        end;
    end;
    amplitudes=mat\B;
 %---------------------------
 %рисование результата разложения
 figure(figure_number);
 plot(1:n,g,'k--');
 hold on;
 for j=1:M
   plot(1:n,amplitudes(j)*phi(j,:),'b-');
   hold on;
 end;
 %общая зависимость
 result=zeros(1,n);
 for i=1:n,
    for j=1:M
    result(i)=result(i)+amplitudes(j)*phi(j,i);
    end;
 end;
plot(1:n,result,'r-');
hold off;
title(strcat('method ',num2str(figure_number)));
else disp('M=0  - there are no gaussians!');
end;