---

import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path


def fourier(signal, fs):
    """Находит частоты, амплитуды и фазы гармоник."""
    n = len(signal)

    coefficients = np.fft.rfft(signal)
    frequencies = np.fft.rfftfreq(n, d=1 / fs)

    # Амплитуды одностороннего спектра
    amplitudes = np.abs(coefficients) / n
    amplitudes[1:] *= 2

    # Частоту Найквиста не удваиваем
    if n % 2 == 0:
        amplitudes[-1] /= 2

    phases = np.angle(coefficients)

    # Выделяем значимые гармоники
    selected = amplitudes > amplitudes.max() * 1e-6

    return (
        frequencies[selected],
        amplitudes[selected],
        phases[selected]
    )


def restore(t, frequencies, amplitudes, phases):
    """Восстанавливает сигнал как сумму гармоник."""
    result = np.zeros_like(t)

    for f, amplitude, phase in zip(
        frequencies, amplitudes, phases
    ):
        result += amplitude * np.cos(
            2 * np.pi * f * t + phase
        )

    return result


def save_figure(fig, filename):
    """Сохраняет график и выводит путь к изображению."""
    fig.tight_layout()
    fig.savefig(filename, dpi=300)
    plt.close(fig)

    print("Сохранено:", Path(filename).resolve())


# ==================================================
# ЗАДАНИЕ 1. СИНУС 500 ГЦ
# ==================================================

f = 500
fs1 = 4 * f          # Частота дискретизации: 2000 Гц
duration1 = 5 / f    # Пять периодов: 0,01 секунды

# Для Фурье берём интервал без повторения конечной точки
t1 = np.arange(round(duration1 * fs1)) / fs1
signal1 = np.sin(2 * np.pi * f * t1)

# Для графика добавляем конечный отсчёт
plot_t1 = np.append(t1, duration1)
plot_signal1 = np.sin(2 * np.pi * f * plot_t1)

# Плотная сетка для отображения исходного синуса
fine_t1 = np.linspace(0, duration1, 2001)
fine_signal1 = np.sin(2 * np.pi * f * fine_t1)

fig1, ax = plt.subplots(figsize=(10, 4))

ax.plot(
    fine_t1 * 1000,
    fine_signal1,
    label="Исходный синус"
)

ax.plot(
    plot_t1 * 1000,
    plot_signal1,
    "o--",
    label=f"Отсчёты: {fs1} Гц"
)

ax.set_title("Задание 1. Синус 500 Гц")
ax.set_xlabel("Время, мс")
ax.set_ylabel("Амплитуда, усл. ед.")
ax.grid()
ax.legend()

save_figure(fig1, "task1_sinus.png")


# ==================================================
# ЗАДАНИЕ 2. АМПЛИТУДНАЯ МОДУЛЯЦИЯ
# ==================================================

fc = 1_000_000       # Частота несущей: 1 МГц
fm = 10_000          # Частота модулирующего сигнала: 10 кГц
A = 1               # Амплитуда несущей
m = 0.5             # Глубина модуляции: 50%
fs2 = 20_000_000     # Частота дискретизации: 20 МГц
duration2 = 3 / fm   # Три периода модулирующего сигнала

t2 = np.arange(round(duration2 * fs2)) / fs2

message = np.sin(2 * np.pi * fm * t2)
carrier = np.cos(2 * np.pi * fc * t2)

envelope = A * (1 + m * message)
signal2 = envelope * carrier

t_us = t2 * 1e6  # Время в микросекундах

fig2, axes = plt.subplots(3, 1, figsize=(11, 9))

# Модулирующий сигнал
axes[0].plot(t_us, message)
axes[0].set_title("Модулирующий сигнал: 10 кГц")

# АМ-сигнал и огибающие
axes[1].plot(
    t_us, signal2,
    linewidth=0.6,
    label="АМ-сигнал"
)
axes[1].plot(t_us, envelope, "r--", label="Огибающие")
axes[1].plot(t_us, -envelope, "r--")

axes[1].set_title(
    "Задание 2. Несущая 1 МГц, глубина модуляции 50%"
)
axes[1].legend()

# Увеличенный фрагмент
fragment = t2 <= 10e-6

axes[2].plot(
    t_us[fragment],
    signal2[fragment],
    "o-",
    markersize=2
)
axes[2].plot(
    t_us[fragment], envelope[fragment], "r--"
)
axes[2].plot(
    t_us[fragment], -envelope[fragment], "r--"
)
axes[2].set_title("Увеличенный фрагмент: первые 10 мкс")

for ax in axes:
    ax.set_xlabel("Время, мкс")
    ax.set_ylabel("Амплитуда, усл. ед.")
    ax.grid()

save_figure(fig2, "task2_modulation.png")


# ==================================================
# ЗАДАНИЕ 3. РАЗЛОЖЕНИЕ ФУРЬЕ
# ==================================================

# Находим гармоники обоих сигналов
frequencies1, amplitudes1, phases1 = fourier(signal1, fs1)
frequencies2, amplitudes2, phases2 = fourier(signal2, fs2)

# Восстанавливаем сигналы суммированием гармоник
restored1 = restore(
    fine_t1, frequencies1, amplitudes1, phases1
)
restored2 = restore(
    t2, frequencies2, amplitudes2, phases2
)

fig3, axes = plt.subplots(2, 2, figsize=(13, 8))

# Восстановление синуса
axes[0, 0].plot(
    plot_t1 * 1000,
    plot_signal1,
    "o",
    label="Исходные отсчёты"
)
axes[0, 0].plot(
    fine_t1 * 1000,
    restored1,
    label="Сумма гармоник"
)
axes[0, 0].set_title("Восстановление синуса")
axes[0, 0].set_xlabel("Время, мс")
axes[0, 0].legend()

# Спектр синуса
axes[0, 1].stem(
    frequencies1, amplitudes1, basefmt="k-"
)
axes[0, 1].set_title("Гармоники синуса")
axes[0, 1].set_xlabel("Частота, Гц")
axes[0, 1].set_xlim(0, fs1 / 2)

# Восстановление АМ-сигнала
axes[1, 0].plot(
    t_us, signal2,
    label="Исходный АМ-сигнал"
)
axes[1, 0].plot(
    t_us, restored2,
    "--",
    label="Сумма гармоник"
)
axes[1, 0].set_title("Восстановление АМ-сигнала")
axes[1, 0].set_xlabel("Время, мкс")
axes[1, 0].set_xlim(0, 10)
axes[1, 0].legend()

# Спектр АМ-сигнала
axes[1, 1].stem(
    frequencies2 / 1000,
    amplitudes2,
    basefmt="k-"
)
axes[1, 1].set_title("Гармоники АМ-сигнала")
axes[1, 1].set_xlabel("Частота, кГц")
axes[1, 1].set_xlim(980, 1020)

for ax in axes.flat:
    ax.set_ylabel("Амплитуда, усл. ед.")
    ax.grid()

fig3.suptitle("Задание 3. Разложение Фурье")

save_figure(fig3, "task3_fourier.png")


# ==================================================
# РЕЗУЛЬТАТЫ В КОНСОЛИ
# ==================================================

print("\nЗадание 1:")
print(f"Частота синуса: {f} Гц")
print(f"Частота дискретизации: {fs1} Гц")
print(f"Шаг дискретизации: {1 / fs1:.6f} с")

print("\nЗадание 2:")
print(f"Частота несущей: {fc} Гц")
print(f"Частота модулирующего сигнала: {fm} Гц")
print(f"Глубина модуляции: {m * 100:.0f}%")
print(f"Минимальная амплитуда: {A * (1 - m)}")
print(f"Максимальная амплитуда: {A * (1 + m)}")

print("\nЗадание 3:")

for name, frequencies, amplitudes, phases in [
    ("Синус", frequencies1, amplitudes1, phases1),
    ("АМ-сигнал", frequencies2, amplitudes2, phases2)
]:
    print("\n" + name)

    for frequency, amplitude, phase in zip(
        frequencies, amplitudes, phases
    ):
        print(
            f"Частота: {frequency:.0f} Гц; "
            f"амплитуда: {amplitude:.3f}; "
            f"фаза: {phase:.3f} рад"
        )

# Проверяем совпадение на исходных отсчётах
error1 = np.max(np.abs(
    signal1 - restore(
        t1, frequencies1, amplitudes1, phases1
    )
))
error2 = np.max(np.abs(signal2 - restored2))

print(f"\nПогрешность восстановления синуса: {error1:.3e}")
print(f"Погрешность восстановления АМ-сигнала: {error2:.3e}")