Преобразование Фурье. Разложение в ряд. Без fft - Matlab - Ответ 8977852

% Период сигнала:
ts=str2double(get(handles.edit1,'String')); 
% Центральная частота, Гц:
f0=str2double(get(handles.edit2,'String')); 
% Девиация частоты, Гц:
dfs=str2double(get(handles.edit3,'String')); 
% Амплитуда сигнала:
A0=str2double(get(handles.edit4,'String')); 
% фаза:
fi0=str2double(get(handles.edit5,'String')); 
% Частота дискретизации
fd=str2double(get(handles.edit6,'String')); 
% Разрешенеи по частоте:
df=str2double(get(handles.edit7,'String')); 
%%
 
% Нижняя частота, Гц:
fn=f0-dfs/2; 
% Верхняя частота, Гц:
fv=f0+dfs/2; 
% Пересчитаем частоты в рад/с:
w0=2*pi*f0; 
wn=2*pi*fn; 
% Вычислим шаг дискретизации:
dt=1/fd; 
%размерность ПФ:
N2=65536;
% Вектор времени:
t=0:dt:ts-dt; 
% ЛЧМ сигнал:
S=A0*exp(-1i*(wn*t+0.5*2*pi*dfs/ts*t.^2+fi0)); 
% Кол-во периодов:
Np=1;
 
for i=1:Np
S(1,(i-1)*length(S)+1:length(S)*i)=S;
end
t=0:dt:Np*ts-dt; 
 
% Временная реализация ЛЧМ сигнала:
% Делаем активными первый оси:
axes(handles.axes1);
% Строим график:
plot(t,real(S))
% Для более наглядного отображения, изменяем пределы по оси У
ylim([-2*A0, 2*A0])
% Включаем сетку:
grid on
% Выводим подписи:
xlabel('Время, с')
ylabel('Амплитуда')
% Включаем зум:
zoom on

% Спектр без fft
axes(handles.axes3); 
% Количество точек:
st=65536;  
% Значения для несущей:
Sn=sin(2*pi*f0*(0:st:2*pi-st));
Cn=cos(2*pi*f0*(0:st:2*pi-st));
 
% выходных переменных:
n=1;
N=65536; 
A=0;
B=0;
% Перебираем частоты:
for F=f0-dfs/2:dfs/N:f0+dfs/2-dfs/N;
    % Находим значения sin и cos для текущих частот:
    Sf=sin(2*pi*F*(0:st:2*pi-st));
    Cf=cos(2*pi*F*(0:st:2*pi-st));
    % Находим А и B
    A(n)=sum(Sn.*Sf);
    B(n)=sum(Cn.*Cf);
    % Добавляем счетчик:
    n=n+1;
end
% Строим график:
plot((f0-dfs/2:dfs/N:f0+dfs/2-dfs/N),sqrt(A.^2+B.^2))
grid on
xlabel('Частота, Гц')
ylabel('АЧХ')
zoom on