Спектральный анализ функции

%Шаг 1. Генерируем шкалу времени от 0 до 1 секунды с шагом в 1/16384 секунды (получится всего 16385 шагов). Массив t из 16385 элементов содержит значения времени.
t = 0:1/16384:1;
 
%Шаг 2. Генерим массивы амплитудами сигналов. В результате получим массивы U1..U5, в каждом из которых 16385 значения сигнала.
U1      = 1             ;
U2      = 2.5  *cos(  t);
U3      = -1.5 *sin(  t);
U4      = 0.7  *cos(2*t);
U5      = -0.2 *sin(3*t);
 
%Шаг 3. Суммируем все сигналы между собой. В результате получим массив x из 16385 значений.
x = U1+U2+U3+U4+U5;
 
%Шаг 4. Для того, чтобы сигнал отображался "красиво", вводим оконную функцию Хемминга. В результате получим массив оконного распределения w из 16385 значений. Операция ".'" нужна для транспонирования оконного массива (он был столбцом, нужно получить строку).
w = hamming(16385).';
 
%Шаг 5. Чтобы было "совсем реально", добавим немного шума. Здесь "size(t)" соответствует "1, 16385", т. е. вызов randn(size(t)) это то же самое, что вызов randn(1, 16385).
 x = x + 0.001*randn(size(t));
 
%Шаг 6. Применим оконную функцию к сигналу x, в результате получим массив z из 16385 значений. Здесь операция ".*" означает "перемножить каждый элемент массива w с соответствующим элементом массива x"
z = w.*x;
 
%Шаг 7. Делаем дискретное преобразование фурье (ДПФ) от сигнала z по 16384 точкам. В результате получим массив комплексных чисел Y (всего 16384 комплексных числа)
Y = fft(z,16384);
 
%Шаг 8. Вычисляем спектральную плотность по мощности, короче говоря - спектр (чего мы и добивались). В результате получим массив Pyy из 16384 значений. Здесь коэффициенты 5000 и 10000 подобраны таким образом, чтобы сместить значение десятичного логагифма в диапазон значений целого беззнакового 16-битного числа 0..65535.
Pyy = 5000 * log10(Y.*conj(Y) * 10000);
 
%Шаг 9. Нарисуем график нашего спектра
%Подготовим ось X - ось частот от 0 до 8191 Гц
f =  (0:8191);
% и нарисуем график. Рисуем не весь массив, а только половину (вторая половина Pyy является зеркальным отображением первой).
plot(f,Pyy(1:8192))
title ('Удельная мощность спектра')
xlabel('Частота (Hz)')
grid

Fs = 1; % Частота дискретизации
T = 1/Fs; % период дискретизации
L = 1000; % длина сигнала
t = (0:L-1)*T; % вектор времени 
x = 1+2.5*cos(t)-1.5*sin(t)+0.7*cos(2*t)-0.2*sin(3*t);
y = x + 0.6*randn(size(t)); % Добавление шума в функцию
subplot(1,2,1)
plot(Fs*t(1:50),y(1:50))
title('Сигнал поврежденный случайным шумом')
xlabel('Время (в миллисекундах)')
 
NFFT = 2^nextpow2(L); % следующая степень числа 2 от длины y
Y = fft(y,NFFT)/L; % fft - одномерное прямое дискретное 
                   % преобразование в ряд Фурье                   
f = Fs/2*linspace(0,1,NFFT/2+1); % linspace - генерация 
                                 % арифметической прогрессии
% Участок одностороннего амплитудного спектра.
subplot(1,2,2)
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Односторонняя амплитуда спектра y(t)')
xlabel('Частота (Гц)')
ylabel('|Y(f)|')

function dxdt = diffsys(t, x)
b=0;
y=0.5;
P=0;
w=2.0;
dxdt(1,1) = x(2); 
dxdt(2,1) = P*cos(w*t)-b*x(2)-x(1)-y*x(1).^3; 
end

clear, clc
 
tk = [0 15]; 
x0 = [0.01 0];
[T X] = ode45('diffsys',tk,x0);
X=X';
for i=1:size(T,1)
    x(i)=X(2,i);
end
L=size(x,2);
NFFT = 2^nextpow2(L); 
Y1 = fft(x,NFFT)/L;                  
f = linspace(0,1,NFFT/2+1);
Y(1)=abs(Y1(1));
for i=2:NFFT/2+1
Y(i)=2*abs(Y1(i));
end
[pks locks]=findpeaks(Y);
stem (f(locks(1:10)), pks(1:10))
set(gca,'XTick',f(locks))
set(gca,'XTickLabel',f(locks)*2*pi)
xlabel('Угловая скорость')
ylabel('|Y(f)|')