%% Моделирование биквадратных фильтров для АЦП с автоматическим расчетом коэффициентов
clear all; close all; clc;

%% Параметры моделирования
Fs = 100000;              % Частота дискретизации [Гц]
T = 0.5;                  % Время моделирования [с]
t = 0:1/Fs:T-1/Fs;      % Временной вектор
N = length(t);          % Количество отсчетов

%% Уровни шума для разных каналов
noise_levels.voltage = 0.2;    % 2% шума для напряжений
noise_levels.current = 0.2;    % 5% шума для токов  
noise_levels.temperature = 0.2; % 1% шума для температур

%% АВТОМАТИЧЕСКИЙ РАСЧЕТ КОЭФФИЦИЕНТОВ ЧЕРЕЗ КЛАСС
fprintf('=== АВТОМАТИЧЕСКИЙ РАСЧЕТ КОЭФФИЦИЕНТОВ ФИЛЬТРОВ ===\n\n');

% 1. Полосовой фильтр 45-55 Гц для напряжений
% [b_bpf, a_bpf, coeffs_bpf] = BiquadFilterDesigner.bpf(20, 10, Fs);
[b_bpf, a_bpf, coeffs_bpf] = BiquadFilterDesigner.lpf(100, Fs);
fprintf('1. Полосовой фильтр 45-55 Гц:\n');
BiquadFilterDesigner.generate_c_code(coeffs_bpf, 'voltage_bpf');

% 2. ФНЧ 100 Гц для токов
[b_lpf_current, a_lpf_current, coeffs_lpf_current] = BiquadFilterDesigner.lpf(100, Fs);
fprintf('2. ФНЧ 100 Гц (токи):\n');
BiquadFilterDesigner.generate_c_code(coeffs_lpf_current, 'current_lpf');

% 3. ФНЧ 10 Гц для температур
[b_lpf_temp, a_lpf_temp, coeffs_lpf_temp] = BiquadFilterDesigner.lpf(10, Fs);
fprintf('3. ФНЧ 10 Гц (температуры):\n');
BiquadFilterDesigner.generate_c_code(coeffs_lpf_temp, 'temperature_lpf');

% Вывод коэффициентов в консоль
fprintf('\n=== РАСЧЕТНЫЕ КОЭФФИЦИЕНТЫ ===\n');
fprintf('Напряжение (BPF 45-55 Гц): b = [%.6f, %.6f, %.6f], a = [1, %.6f, %.6f]\n', ...
        b_bpf, a_bpf(2), a_bpf(3));
fprintf('Ток (LPF 100 Гц): b = [%.6f, %.6f, %.6f], a = [1, %.6f, %.6f]\n', ...
        b_lpf_current, a_lpf_current(2), a_lpf_current(3));
fprintf('Температура (LPF 10 Гц): b = [%.6f, %.6f, %.6f], a = [1, %.6f, %.6f]\n\n', ...
        b_lpf_temp, a_lpf_temp(2), a_lpf_temp(3));

%% Генерация тестовых сигналов

% 1. Сетевое напряжение 50 Гц + гармоники
f_net = 50;
voltage_clean = 1.0 * sin(2*pi*f_net*t) + ...      % основная 50 Гц
                0.1 * sin(2*pi*3*f_net*t) + ...    % 3-я гармоника 150 Гц
                0.05 * sin(2*pi*5*f_net*t);        % 5-я гармоника 250 Гц

% 2. Ток с модуляцией
f_current = 50;
current_clean = 1 * sin(2*pi*f_current*t); %.* ...
                %(1 + 0.2 * sin(2*pi*2*t));         % амплитудная модуляция 2 Гц

% 3. Температура (медленно меняющийся сигнал)
temperature_clean = 25 + 2 * sin(2*pi*0.1*t) + ... % медленные колебания 0.1 Гц
                    0.5 * sin(2*pi*1*t);           % быстрые колебания 1 Гц

%% Добавление шума
voltage_noisy = voltage_clean + noise_levels.voltage * randn(size(t));
current_noisy = current_clean + noise_levels.current * randn(size(t));  
temperature_noisy = temperature_clean + noise_levels.temperature * randn(size(t));

%% Фильтрация сигналов
voltage_filtered = filter(b_bpf, a_bpf, voltage_noisy);
current_filtered = filter(b_lpf_current, a_lpf_current, current_noisy);
temperature_filtered = filter(b_lpf_temp, a_lpf_temp, temperature_noisy);

%% НОРМАЛИЗАЦИЯ УСИЛЕНИЯ (важно для правильного SNR)
% Получаем АЧХ для нормализации
[h_bpf, f_bpf] = freqz(b_bpf, a_bpf, 1024, Fs);
[h_lpf_curr, f_lpf_curr] = freqz(b_lpf_current, a_lpf_current, 1024, Fs);
[h_lpf_temp, f_lpf_temp] = freqz(b_lpf_temp, a_lpf_temp, 1024, Fs);

% Нормализация полосового фильтра на центральной частоте
[~, idx_50hz] = min(abs(f_bpf - 50));
gain_bpf = abs(h_bpf(idx_50hz));
if gain_bpf > 0
    voltage_filtered = voltage_filtered / gain_bpf;
end

% Нормализация ФНЧ на постоянной составляющей
gain_lpf_current = sum(b_lpf_current) / (1 + sum(a_lpf_current(2:end)));
if gain_lpf_current > 0
    current_filtered = current_filtered / gain_lpf_current;
end

gain_lpf_temp = sum(b_lpf_temp) / (1 + sum(a_lpf_temp(2:end)));
if gain_lpf_temp > 0
    temperature_filtered = temperature_filtered / gain_lpf_temp;
end

%% ОКНО 1: НАПРЯЖЕНИЕ
figure('Name', 'Анализ напряжения');

% Временные характеристики
subplot(2,3,1);
plot(t, voltage_noisy, 'b', 'LineWidth', 1); hold on;
plot(t, voltage_filtered, 'r', 'LineWidth', 2);
plot(t, voltage_clean, 'g--', 'LineWidth', 1);
title('Напряжение: временная область');
legend('С шумом', 'Фильтрованный', 'Идеальный', 'Location', 'best');
xlabel('Время [с]'); ylabel('Амплитуда');
grid on;

% АЧХ фильтра
subplot(2,3,2);
plot(f_bpf, 20*log10(abs(h_bpf)), 'LineWidth', 2);
title('АЧХ: Полосовой фильтр 45-55 Гц');
xlabel('Частота [Гц]'); ylabel('Усиление [дБ]');
grid on; xlim([0, 200]);

% Спектр фильтрованного сигнала
subplot(2,3,3);
[P_voltage, f_voltage] = pwelch(voltage_filtered, [], [], 1024, Fs);
plot(f_voltage, 10*log10(P_voltage), 'LineWidth', 2);
title('Спектр фильтрованного напряжения');
xlabel('Частота [Гц]'); ylabel('Мощность [дБ]');
grid on; xlim([0, 200]);

% Групповая задержка
subplot(2,3,4);
[gd_bpf, f_gd] = grpdelay(b_bpf, a_bpf, 1024, Fs);
plot(f_gd, gd_bpf/Fs*1000, 'LineWidth', 2);
title('Групповая задержка фильтра');
xlabel('Частота [Гц]'); ylabel('Задержка [мс]');
grid on; xlim([0, 100]);

% Детальный вид (последние 0.1 секунды)
idx_end_voltage = max(1, length(t) - 0.1*Fs + 1):length(t); % последние 100 мс
subplot(2,3,5);
plot(t(idx_end_voltage), voltage_noisy(idx_end_voltage), 'b', 'LineWidth', 2); hold on;
plot(t(idx_end_voltage), voltage_filtered(idx_end_voltage), 'r', 'LineWidth', 2);
plot(t(idx_end_voltage), voltage_clean(idx_end_voltage), 'g--', 'LineWidth', 2);
title('Напряжение: УВЕЛИЧЕННЫЙ ВИД (900-1000 мс)');
legend('С шумом', 'Фильтрованный', 'Идеальный', 'Location', 'best');
xlabel('Время [с]'); ylabel('Амплитуда');
grid on;
xlim([t(idx_end_voltage(1)) t(idx_end_voltage(end))]);

% Статистика
subplot(2,3,6);
snr_voltage_in = snr(voltage_clean, voltage_noisy - voltage_clean);
snr_voltage_out = snr(voltage_clean, voltage_filtered - voltage_clean);
improvement_voltage = snr_voltage_out - snr_voltage_in;

text(0.1, 0.8, sprintf('SNR вход: %.1f дБ', snr_voltage_in), 'FontSize', 12);
text(0.1, 0.6, sprintf('SNR выход: %.1f дБ', snr_voltage_out), 'FontSize', 12);
text(0.1, 0.4, sprintf('Улучшение: %.1f дБ', improvement_voltage), 'FontSize', 12);
text(0.1, 0.2, sprintf('Задержка 50 Гц: %.1f мс', gd_bpf(idx_50hz)/Fs*1000), 'FontSize', 12);
title('Статистика фильтрации напряжения');
axis off;

%% ОКНО 2: ТОК
figure('Name', 'Анализ тока');

% Временные характеристики
subplot(2,3,1);
plot(t, current_noisy, 'b', 'LineWidth', 1); hold on;
plot(t, current_filtered, 'r', 'LineWidth', 2);
plot(t, current_clean, 'g--', 'LineWidth', 1);
title('Ток: временная область');
legend('С шумом', 'Фильтрованный', 'Идеальный', 'Location', 'best');
xlabel('Время [с]'); ylabel('Амплитуда');
grid on;

% АЧХ фильтра
subplot(2,3,2);
plot(f_lpf_curr, 20*log10(abs(h_lpf_curr)), 'LineWidth', 2);
title('АЧХ: ФНЧ 100 Гц');
xlabel('Частота [Гц]'); ylabel('Усиление [дБ]');
grid on; xlim([0, 200]);

% Спектр фильтрованного сигнала
subplot(2,3,3);
[P_current, f_current] = pwelch(current_filtered, [], [], 1024, Fs);
plot(f_current, 10*log10(P_current), 'LineWidth', 2);
title('Спектр фильтрованного тока');
xlabel('Частота [Гц]'); ylabel('Мощность [дБ]');
grid on; xlim([0, 200]);

% Групповая задержка
subplot(2,3,4);
[gd_lpf_curr, f_gd] = grpdelay(b_lpf_current, a_lpf_current, 1024, Fs);
plot(f_gd, gd_lpf_curr/Fs*1000, 'LineWidth', 2);
title('Групповая задержка фильтра');
xlabel('Частота [Гц]'); ylabel('Задержка [мс]');
grid on; xlim([0, 100]);

% Детальный вид (последние 0.1 секунды)
idx_end_current = max(1, length(t) - 0.1*Fs + 1):length(t); % последние 100 мс
subplot(2,3,5);
plot(t(idx_end_current), current_noisy(idx_end_current), 'b', 'LineWidth', 2); hold on;
plot(t(idx_end_current), current_filtered(idx_end_current), 'r', 'LineWidth', 2);
plot(t(idx_end_current), current_clean(idx_end_current), 'g--', 'LineWidth', 2);
title('Ток: УВЕЛИЧЕННЫЙ ВИД (900-1000 мс)');
legend('С шумом', 'Фильтрованный', 'Идеальный', 'Location', 'best');
xlabel('Время [с]'); ylabel('Амплитуда');
grid on;
xlim([t(idx_end_current(1)) t(idx_end_current(end))]);

% Статистика
subplot(2,3,6);
snr_current_in = snr(current_clean, current_noisy - current_clean);
snr_current_out = snr(current_clean, current_filtered - current_clean);
improvement_current = snr_current_out - snr_current_in;

[~, idx_50hz_curr] = min(abs(f_gd - 50));
text(0.1, 0.8, sprintf('SNR вход: %.1f дБ', snr_current_in), 'FontSize', 12);
text(0.1, 0.6, sprintf('SNR выход: %.1f дБ', snr_current_out), 'FontSize', 12);
text(0.1, 0.4, sprintf('Улучшение: %.1f дБ', improvement_current), 'FontSize', 12);
text(0.1, 0.2, sprintf('Задержка 50 Гц: %.1f мс', gd_lpf_curr(idx_50hz_curr)/Fs*1000), 'FontSize', 12);
title('Статистика фильтрации тока');
axis off;

%% ОКНО 3: ТЕМПЕРАТУРА
figure('Name', 'Анализ температуры');

% Временные характеристики
subplot(2,3,1);
plot(t, temperature_noisy, 'b', 'LineWidth', 1); hold on;
plot(t, temperature_filtered, 'r', 'LineWidth', 2);
plot(t, temperature_clean, 'g--', 'LineWidth', 1);
title('Температура: временная область');
legend('С шумом', 'Фильтрованный', 'Идеальный', 'Location', 'best');
xlabel('Время [с]'); ylabel('Температура [°C]');
grid on;

% АЧХ фильтра
subplot(2,3,2);
plot(f_lpf_temp, 20*log10(abs(h_lpf_temp)), 'LineWidth', 2);
title('АЧХ: ФНЧ 10 Гц');
xlabel('Частота [Гц]'); ylabel('Усиление [дБ]');
grid on; xlim([0, 20]);

% Спектр фильтрованного сигнала
subplot(2,3,3);
[P_temp, f_temp] = pwelch(temperature_filtered, [], [], 1024, Fs);
plot(f_temp, 10*log10(P_temp), 'LineWidth', 2);
title('Спектр фильтрованной температуры');
xlabel('Частота [Гц]'); ylabel('Мощность [дБ]');
grid on; xlim([0, 10]);

% Групповая задержка
subplot(2,3,4);
[gd_lpf_temp, f_gd_temp] = grpdelay(b_lpf_temp, a_lpf_temp, 1024, Fs);
plot(f_gd_temp, gd_lpf_temp/Fs*1000, 'LineWidth', 2);
title('Групповая задержка фильтра');
xlabel('Частота [Гц]'); ylabel('Задержка [мс]');
grid on; xlim([0, 20]);

% Детальный вид (последне 0.1 секунды)
idx_end_temp = max(1, length(t) - 0.1*Fs + 1):length(t); % последние 100 мс
subplot(2,3,5);
plot(t(idx_end_temp), temperature_noisy(idx_end_temp), 'b', 'LineWidth', 2); hold on;
plot(t(idx_end_temp), temperature_filtered(idx_end_temp), 'r', 'LineWidth', 2);
plot(t(idx_end_temp), temperature_clean(idx_end_temp), 'g--', 'LineWidth', 2);
title('Температура: УВЕЛИЧЕННЫЙ ВИД (900-1000 мс)');
legend('С шумом', 'Фильтрованный', 'Идеальный', 'Location', 'best');
xlabel('Время [с]'); ylabel('Температура [°C]');
grid on;
xlim([t(idx_end_temp(1)) t(idx_end_temp(end))]);

% Статистика
subplot(2,3,6);
snr_temp_in = snr(temperature_clean, temperature_noisy - temperature_clean);
snr_temp_out = snr(temperature_clean, temperature_filtered - temperature_clean);
improvement_temp = snr_temp_out - snr_temp_in;

[~, idx_1hz] = min(abs(f_gd_temp - 1));
text(0.1, 0.8, sprintf('SNR вход: %.1f дБ', snr_temp_in), 'FontSize', 12);
text(0.1, 0.6, sprintf('SNR выход: %.1f дБ', snr_temp_out), 'FontSize', 12);
text(0.1, 0.4, sprintf('Улучшение: %.1f дБ', improvement_temp), 'FontSize', 12);
text(0.1, 0.2, sprintf('Задержка 1 Гц: %.1f мс', gd_lpf_temp(idx_1hz)/Fs*1000), 'FontSize', 12);
title('Статистика фильтрации температуры');
axis off;

%% Вывод результатов в командное окно
fprintf('\n=== ИТОГИ ФИЛЬТРАЦИИ С АВТОРАСЧЕТОМ КОЭФФИЦИЕНТОВ ===\n\n');
fprintf('НАПРЯЖЕНИЕ (полосовой 45-55 Гц):\n');
fprintf('  SNR вход: %.1f дБ, выход: %.1f дБ, улучшение: %.1f дБ\n', ...
        snr_voltage_in, snr_voltage_out, improvement_voltage);
fprintf('  Задержка на 50 Гц: %.1f мс\n\n', gd_bpf(idx_50hz)/Fs*1000);

fprintf('ТОК (ФНЧ 100 Гц):\n');
fprintf('  SNR вход: %.1f дБ, выход: %.1f дБ, улучшение: %.1f дБ\n', ...
        snr_current_in, snr_current_out, improvement_current);
fprintf('  Задержка на 50 Гц: %.1f мс\n\n', gd_lpf_curr(idx_50hz_curr)/Fs*1000);

fprintf('ТЕМПЕРАТУРА (ФНЧ 10 Гц):\n');
fprintf('  SNR вход: %.1f дБ, выход: %.1f дБ, улучшение: %.1f дБ\n', ...
        snr_temp_in, snr_temp_out, improvement_temp);
fprintf('  Задержка на 1 Гц: %.1f мс\n\n', gd_lpf_temp(idx_1hz)/Fs*1000);