---

// ====================================================================
// Лабораторная работа: ДПФ и ОДПФ в Scilab
// ====================================================================
clear;
clc;
clf;

// 1. Исходные данные с доски
N  = 8;             // Количество отсчетов
fs = 8000;          // Частота дискретизации, Гц
f1 = 1000;          // Частота 1-й гармоники, Гц
f2 = 2000;          // Частота 2-й гармоники, Гц

// Амплитуды и фазы (при необходимости скорректируйте под свой вариант)
A1 = 1.0;
A2 = 1.0;
phi1 = 0.0;
phi2 = 0.0;

ts = 1 / fs;        // Период дискретизации, с
n_idx = 0:(N-1);    // Индексы отсчетов n = 0..7
t = n_idx * ts;     // Вектор времени t = 0:ts:(N-1)*ts

// 2. Формирование дискретного сигнала x(n)
x = A1 * sin(2 * %pi * f1 * t + phi1) + A2 * sin(2 * %pi * f2 * t + phi2);

mprintf("------------------- ОТСЧЕТЫ СИГНАЛА x(n) -------------------\n");
for n = 0:(N-1)
    mprintf("n = %d | t = %.6f с | x(n) = %8.4f\n", n, t(n+1), x(n+1));
end

// ====================================================================
// 3. Задание 1: Вычисление ДПФ вручную по формуле и сравнение с fft
//    X(m) = sum_{n=0}^{N-1} x(n) * [cos(2*pi*n*m/N) - %i * sin(2*pi*n*m/N)]
// ====================================================================
X_manual = zeros(1, N);

for m = 0:(N-1)
    sum_val = 0;
    for n = 0:(N-1)
        angle = 2 * %pi * n * m / N;
        sum_val = sum_val + x(n+1) * (cos(angle) - %i * sin(angle));
    end
    X_manual(m+1) = sum_val;
end

// Встроенное ДПФ (в Scilab fft(x, -1) реализует прямое преобразование)
X_fft = fft(x, -1);

mprintf("\n------------------- ЗАДАНИЕ 1: РЕЗУЛЬТАТЫ ДПФ -------------------\n");
mprintf("m | Частота, Гц |    X_manual (по формуле)    |        X_fft (встроенное)\n");
mprintf("--+-------------+-----------------------------+--------------------------\n");
for m = 0:(N-1)
    fm = m * fs / N;
    mprintf("%d |   %5.0f     |  %7.3f %+7.3fj        |  %7.3f %+7.3fj\n", ..
            m, fm, real(X_manual(m+1)), imag(X_manual(m+1)), ..
                   real(X_fft(m+1)), imag(X_fft(m+1)));
end

err_dft = max(abs(X_manual - X_fft));
mprintf("\nМаксимальное расхождение между формулой и fft: %.2e\n", err_dft);

// ====================================================================
// 4. Задание 2: Вычисление ОДПФ вручную по формуле и сравнение с ifft
//    x(n) = (1/N) * sum_{m=0}^{N-1} X(m) * [cos(2*pi*n*m/N) + %i * sin(2*pi*n*m/N)]
// ====================================================================
x_rec_manual = zeros(1, N);

for n = 0:(N-1)
    sum_val = 0;
    for m = 0:(N-1)
        angle = 2 * %pi * n * m / N;
        sum_val = sum_val + X_manual(m+1) * (cos(angle) + %i * sin(angle));
    end
    x_rec_manual(n+1) = sum_val / N;
end

// Встроенное обратное ДПФ (fft(..., 1) или ifft(...))
if isdef("ifft") then
    x_ifft = ifft(X_manual);
else
    x_ifft = fft(X_manual, 1);
end

mprintf("\n------------------- ЗАДАНИЕ 2: РЕЗУЛЬТАТЫ ОДПФ -------------------\n");
mprintf("n | Исходный x(n) | x_rec (по формуле) | x_ifft (встроенное)\n");
mprintf("--+---------------+--------------------+--------------------\n");
for n = 0:(N-1)
    mprintf("%d |    %8.4f   |       %8.4f     |       %8.4f\n", ..
            n, x(n+1), real(x_rec_manual(n+1)), real(x_ifft(n+1)));
end

err_idft = max(abs(x - real(x_rec_manual)));
mprintf("\nМаксимальная ошибка восстановления сигнала через ОДПФ: %.2e\n", err_idft);

// ====================================================================
// 5. Построение графиков
// ====================================================================
// 1) Исходный дискретный сигнал
subplot(2, 2, 1);
plot2d3(n_idx, x, -9);
plot(n_idx, x, "r.");
xtitle("Исходный дискретный сигнал x(n)", "Номер отсчета n", "Амплитуда");
xgrid();

// 2) Амплитудный спектр |X(m)|
subplot(2, 2, 2);
f_axis = n_idx * fs / N;
plot2d3(f_axis, abs(X_manual), -9);
plot(f_axis, abs(X_manual), "b.");
xtitle("Амплитудный спектр |X(m)|", "Частота, Гц", "|X(m)|");
xgrid();

// 3) Фазовый спектр arg(X(m)) (для значимых спектральных линий)
subplot(2, 2, 3);
phase_spectrum = atan(imag(X_manual), real(X_manual));
// Обнуление фазы там, где модуль пренебрежимо мал
for k = 1:N
    if abs(X_manual(k)) < 1e-6 then
        phase_spectrum(k) = 0;
    end
end
plot2d3(f_axis, phase_spectrum, -9);
plot(f_axis, phase_spectrum, "g.");
xtitle("Фазовый спектр arg(X(m))", "Частота, Гц", "Фаза, рад");
xgrid();

// 4) Сравнение исходного и восстановленного сигнала
subplot(2, 2, 4);
plot(n_idx, x, "b-o");
plot(n_idx, real(x_rec_manual), "r--x");
legend(["Исходный сигнал x(n)", "Восстановленный ОДПФ"]);
xtitle("Сравнение x(n) и восстановленного сигнала", "Номер отсчета n", "Амплитуда");
xgrid();