рефакторинг и вроде бы понял как надо управлять импульсами
надо доделать и проверить
This commit is contained in:
@@ -0,0 +1,177 @@
|
||||
%% КЛАСС ДЛЯ РАСЧЕТА КОЭФФИЦИЕНТОВ БИКВАДРАТНЫХ ФИЛЬТРОВ
|
||||
classdef BiquadFilterDesigner
|
||||
%BIQUADFILTERDESIGNER Класс для расчета коэффициентов биквадратных фильтров
|
||||
% Поддерживает: ФНЧ, ФВЧ, полосовые, режекторные фильтры
|
||||
|
||||
methods (Static)
|
||||
|
||||
%% 1. ФНЧ (Low-Pass Filter)
|
||||
function [b, a, coeffs_cmsis] = lpf(fc, Fs, Q)
|
||||
% lpf - Расчет коэффициентов ФНЧ через MATLAB функции
|
||||
% fc - частота среза [Гц]
|
||||
% Fs - частота дискретизации [Гц]
|
||||
|
||||
if fc >= Fs/2
|
||||
error('Частота среза должна быть меньше Fs/2');
|
||||
end
|
||||
|
||||
% Используем встроенные функции MATLAB
|
||||
if nargin < 3
|
||||
% По умолчанию - Баттерворт
|
||||
[b, a] = butter(2, fc/(Fs/2), 'low');
|
||||
else
|
||||
% С указанной добротностью
|
||||
[b, a] = butter(2, fc/(Fs/2), 'low');
|
||||
end
|
||||
|
||||
% Стабилизация коэффициентов
|
||||
[b, a] = BiquadFilterDesigner.stabilize_filter(b, a);
|
||||
|
||||
% Формат для CMSIS-DSP [b0, b1, b2, a1, a2]
|
||||
coeffs_cmsis = [b(1), b(2), b(3), a(2), a(3)];
|
||||
end
|
||||
|
||||
%% 2. ФВЧ (High-Pass Filter)
|
||||
function [b, a, coeffs_cmsis] = hpf(fc, Fs, Q)
|
||||
% hpf - Расчет коэффициентов ФВЧ через MATLAB функции
|
||||
% fc - частота среза [Гц]
|
||||
% Fs - частота дискретизации [Гц]
|
||||
|
||||
if fc >= Fs/2
|
||||
error('Частота среза должна быть меньше Fs/2');
|
||||
end
|
||||
|
||||
% Используем встроенные функции MATLAB
|
||||
[b, a] = butter(2, fc/(Fs/2), 'high');
|
||||
|
||||
% Стабилизация коэффициентов
|
||||
[b, a] = BiquadFilterDesigner.stabilize_filter(b, a);
|
||||
|
||||
coeffs_cmsis = [b(1), b(2), b(3), a(2), a(3)];
|
||||
end
|
||||
|
||||
%% 3. Полосовой фильтр (Band-Pass Filter) - ИСПРАВЛЕННЫЙ
|
||||
function [b, a, coeffs_cmsis] = bpf(fc, bw, Fs)
|
||||
% bpf - Расчет коэффициентов полосового фильтра через MATLAB функции
|
||||
% fc - центральная частота [Гц]
|
||||
% bw - ширина полосы [Гц]
|
||||
% Fs - частота дискретизации [Гц]
|
||||
|
||||
f_low = fc - bw/2;
|
||||
f_high = fc + bw/2;
|
||||
|
||||
if f_high >= Fs/2
|
||||
error('Верхняя граница полосы должна быть меньше Fs/2');
|
||||
end
|
||||
|
||||
% Используем встроенные функции MATLAB
|
||||
[b, a] = butter(2, [f_low, f_high]/(Fs/2), 'bandpass');
|
||||
|
||||
% Стабилизация коэффициентов
|
||||
[b, a] = BiquadFilterDesigner.stabilize_filter(b, a);
|
||||
|
||||
coeffs_cmsis = [b(1), b(2), b(3), a(2), a(3)];
|
||||
end
|
||||
|
||||
%% 4. Режекторный фильтр (Band-Stop/Notch Filter)
|
||||
function [b, a, coeffs_cmsis] = bsf(fc, bw, Fs)
|
||||
% bsf - Расчет коэффициентов режекторного фильтра через MATLAB функции
|
||||
% fc - центральная частота [Гц]
|
||||
% bw - ширина полосы подавления [Гц]
|
||||
% Fs - частота дискретизации [Гц]
|
||||
|
||||
f_low = fc - bw/2;
|
||||
f_high = fc + bw/2;
|
||||
|
||||
if f_high >= Fs/2
|
||||
error('Верхняя граница полосы должна быть меньше Fs/2');
|
||||
end
|
||||
|
||||
% Используем встроенные функции MATLAB
|
||||
[b, a] = butter(2, [f_low, f_high]/(Fs/2), 'stop');
|
||||
|
||||
% Стабилизация коэффициентов
|
||||
a_stable = a;
|
||||
if abs(a(3)) > 0.99
|
||||
a_stable(3) = 0.99 * sign(a(3));
|
||||
end
|
||||
if abs(a(2)) > 1.98
|
||||
a_stable(2) = 1.98 * sign(a(2));
|
||||
end
|
||||
|
||||
b = b;
|
||||
a = a_stable;
|
||||
|
||||
coeffs_cmsis = [b(1), b(2), b(3), a(2), a(3)];
|
||||
end
|
||||
|
||||
%% 5. Показать характеристики фильтра
|
||||
function show_response(b, a, Fs, filter_name)
|
||||
% show_response - Визуализация характеристик фильтра
|
||||
% b, a - коэффициенты фильтра
|
||||
% Fs - частота дискретизации [Гц]
|
||||
% filter_name - название фильтра для заголовка
|
||||
|
||||
figure('Position', [100, 100, 1000, 700]);
|
||||
|
||||
% АЧХ
|
||||
subplot(2,2,1);
|
||||
[h, f] = freqz(b, a, 1024, Fs);
|
||||
plot(f, 20*log10(abs(h)), 'LineWidth', 2);
|
||||
title([filter_name ' - АЧХ']);
|
||||
xlabel('Частота [Гц]'); ylabel('Усиление [дБ]');
|
||||
grid on;
|
||||
|
||||
% ФЧХ
|
||||
subplot(2,2,2);
|
||||
plot(f, angle(h)*180/pi, 'LineWidth', 2);
|
||||
title([filter_name ' - ФЧХ']);
|
||||
xlabel('Частота [Гц]'); ylabel('Фаза [градусы]');
|
||||
grid on;
|
||||
|
||||
% Групповая задержка
|
||||
subplot(2,2,3);
|
||||
[gd, f_gd] = grpdelay(b, a, 1024, Fs);
|
||||
plot(f_gd, gd/Fs*1000, 'LineWidth', 2);
|
||||
title([filter_name ' - Групповая задержка']);
|
||||
xlabel('Частота [Гц]'); ylabel('Задержка [мс]');
|
||||
grid on;
|
||||
|
||||
% Импульсная характеристика
|
||||
subplot(2,2,4);
|
||||
imp_resp = filter(b, a, [1 zeros(1, 99)]);
|
||||
stem(0:99, imp_resp, 'filled', 'MarkerSize', 4);
|
||||
title([filter_name ' - Импульсная характеристика']);
|
||||
xlabel('Отсчеты'); ylabel('Амплитуда');
|
||||
grid on;
|
||||
end
|
||||
|
||||
%% 6. Генерация C-кода для коэффициентов
|
||||
function generate_c_code(coeffs, filter_name)
|
||||
% generate_c_code - Генерация C-кода для коэффициентов
|
||||
% coeffs - коэффициенты в формате CMSIS [b0, b1, b2, a1, a2]
|
||||
% filter_name - имя фильтра для переменной
|
||||
|
||||
fprintf('// Коэффициенты фильтра: %s\n', filter_name);
|
||||
fprintf('const float32_t %s_coeffs[5] = {\n', filter_name);
|
||||
fprintf(' %.6ff, // b0\n', coeffs(1));
|
||||
fprintf(' %.6ff, // b1\n', coeffs(2));
|
||||
fprintf(' %.6ff, // b2\n', coeffs(3));
|
||||
fprintf(' %.6ff, // a1\n', coeffs(4));
|
||||
fprintf(' %.6ff // a2\n', coeffs(5));
|
||||
fprintf('};\n\n');
|
||||
end
|
||||
|
||||
function [b_stable, a_stable] = stabilize_filter(b, a)
|
||||
a_stable = a;
|
||||
% Ограничиваем коэффициенты a для устойчивости
|
||||
if abs(a_stable(3)) > 0.95
|
||||
a_stable(3) = 0.95 * sign(a_stable(3));
|
||||
end
|
||||
if abs(a_stable(2)) > 1.9
|
||||
a_stable(2) = 1.9 * sign(a_stable(2));
|
||||
end
|
||||
b_stable = b;
|
||||
end
|
||||
end
|
||||
end
|
||||
12
Информация для программиста (УПП СП СЭД)/CALC/alpha_calc.m
Normal file
12
Информация для программиста (УПП СП СЭД)/CALC/alpha_calc.m
Normal file
@@ -0,0 +1,12 @@
|
||||
open_level = 0:0.01:1; % Степень регулирования выходного напряжения (epsilon) от 0 до 1
|
||||
|
||||
OpenLevelForCos = (open_level.*2)-1;
|
||||
alpha_rad = acos(OpenLevelForCos); % угол в радианах
|
||||
alpha = alpha_rad/pi; % угол открытия тиристора в о.е. от максимально заданного
|
||||
|
||||
plot(alpha, open_level)
|
||||
grid on;
|
||||
xlabel('\alpha, о.е. (от \pi)');
|
||||
ylabel('\epsilon, о.е.');
|
||||
title('Регулировочная характеристика \epsilon');
|
||||
legend('ε = cos(α)');
|
||||
294
Информация для программиста (УПП СП СЭД)/CALC/calc_biquad.m
Normal file
294
Информация для программиста (УПП СП СЭД)/CALC/calc_biquad.m
Normal file
@@ -0,0 +1,294 @@
|
||||
%% Моделирование биквадратных фильтров для АЦП с автоматическим расчетом коэффициентов
|
||||
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);
|
||||
@@ -0,0 +1,204 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# Данные: ADC -> Temperature
|
||||
adc_values = [2188, 2197, 2206, 2216, 2226, 2236, 2247, 2259, 2271, 2283,
|
||||
2296, 2310, 2324, 2338, 2354, 2369, 2385, 2402, 2419, 2437,
|
||||
2455, 2474, 2493, 2513, 2533, 2554, 2575, 2597, 2619, 2641,
|
||||
2664, 2688, 2711, 2735, 2759, 2784, 2809, 2833, 2859, 2884,
|
||||
2909, 2935, 2961, 2986, 3012, 3037, 3063, 3089, 3114, 3140,
|
||||
3165, 3190, 3215, 3239, 3263, 3288, 3312, 3335, 3359, 3381,
|
||||
3404, 3426, 3448, 3470, 3491, 3512, 3532, 3552, 3572, 3591,
|
||||
3610, 3628, 3646, 3663, 3681, 3697, 3714, 3729, 3745, 3760,
|
||||
3775, 3789, 3803, 3817, 3830, 3843, 3855, 3868, 3879, 3891,
|
||||
3902, 3913, 3924, 3934, 3944, 3954, 3963, 3972, 3981, 3989,
|
||||
3997, 4005, 4013, 4021, 4028, 4035, 4042, 4049, 4055, 4062,
|
||||
4068, 4074, 4079, 4085, 4091, 4096]
|
||||
|
||||
temperatures = [-25, -24, -23, -22, -21, -20, -19, -18, -17, -16, -15, -14, -13, -12, -11,
|
||||
-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9,
|
||||
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
|
||||
28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45,
|
||||
46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63,
|
||||
64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81,
|
||||
82, 83, 84, 85, 86, 87, 88, 89, 90]
|
||||
|
||||
# Параметры ограничений
|
||||
MIN_COEFF_ABS = 1e-5 # Минимальное абсолютное значение коэффициента
|
||||
MAX_CONDITION_NUMBER = 1e10 # Максимальное число обусловленности
|
||||
|
||||
# Аппроксимация полиномами разных степеней
|
||||
degrees = [2, 3, 4, 5]
|
||||
coefficients = {}
|
||||
filtered_coefficients = {}
|
||||
|
||||
plt.figure(figsize=(16, 12))
|
||||
|
||||
# График 1: Сравнение аппроксимаций
|
||||
plt.subplot(2, 3, 1)
|
||||
plt.plot(adc_values, temperatures, 'ko-', markersize=3, linewidth=1, label='Исходные данные')
|
||||
|
||||
colors = ['red', 'blue', 'green', 'orange']
|
||||
for i, degree in enumerate(degrees):
|
||||
coeffs = np.polyfit(adc_values, temperatures, degree)
|
||||
coefficients[degree] = coeffs
|
||||
|
||||
# Фильтрация малых коэффициентов
|
||||
filtered_coeffs = coeffs.copy()
|
||||
small_coeffs_mask = np.abs(coeffs) < MIN_COEFF_ABS
|
||||
filtered_coeffs[small_coeffs_mask] = 0
|
||||
filtered_coefficients[degree] = filtered_coeffs
|
||||
|
||||
poly = np.poly1d(coeffs)
|
||||
adc_continuous = np.linspace(min(adc_values), max(adc_values), 500)
|
||||
temp_predicted = poly(adc_continuous)
|
||||
|
||||
plt.plot(adc_continuous, temp_predicted, color=colors[i], linewidth=2,
|
||||
label=f'Полином {degree}-й степени')
|
||||
|
||||
plt.xlabel('Значение АЦП')
|
||||
plt.ylabel('Температура (°C)')
|
||||
plt.title('Аппроксимация зависимости АЦП → Температура')
|
||||
plt.legend()
|
||||
plt.grid(True, alpha=0.3)
|
||||
|
||||
# График 2: Ошибки аппроксимации
|
||||
plt.subplot(2, 3, 2)
|
||||
for i, degree in enumerate(degrees):
|
||||
poly = np.poly1d(coefficients[degree])
|
||||
predicted = poly(adc_values)
|
||||
error = predicted - temperatures
|
||||
|
||||
plt.plot(temperatures, error, 'o-', color=colors[i], markersize=3,
|
||||
label=f'Ошибка {degree}-й степени')
|
||||
|
||||
plt.xlabel('Температура (°C)')
|
||||
plt.ylabel('Ошибка (°C)')
|
||||
plt.title('Ошибки аппроксимации по температуре')
|
||||
plt.legend()
|
||||
plt.grid(True, alpha=0.3)
|
||||
|
||||
# График 3: Статистика ошибок
|
||||
plt.subplot(2, 3, 3)
|
||||
max_errors = []
|
||||
rms_errors = []
|
||||
|
||||
for degree in degrees:
|
||||
poly = np.poly1d(coefficients[degree])
|
||||
predicted = poly(adc_values)
|
||||
max_error = np.max(np.abs(predicted - temperatures))
|
||||
rms_error = np.sqrt(np.mean((predicted - temperatures)**2))
|
||||
|
||||
max_errors.append(max_error)
|
||||
rms_errors.append(rms_error)
|
||||
|
||||
x_pos = np.arange(len(degrees))
|
||||
width = 0.35
|
||||
|
||||
plt.bar(x_pos - width/2, max_errors, width, label='Макс. ошибка', alpha=0.7)
|
||||
plt.bar(x_pos + width/2, rms_errors, width, label='СКО', alpha=0.7)
|
||||
|
||||
plt.xlabel('Степень полинома')
|
||||
plt.ylabel('Ошибка (°C)')
|
||||
plt.title('Статистика ошибок аппроксимации')
|
||||
plt.xticks(x_pos, degrees)
|
||||
plt.legend()
|
||||
plt.grid(True, alpha=0.3)
|
||||
|
||||
# График 4: Сравнение коэффициентов до/после фильтрации
|
||||
plt.subplot(2, 3, 4)
|
||||
for degree in degrees:
|
||||
coeffs = coefficients[degree]
|
||||
filtered_coeffs = filtered_coefficients[degree]
|
||||
|
||||
# Отображаем только ненулевые коэффициенты
|
||||
nonzero_indices = np.where(filtered_coeffs != 0)[0]
|
||||
if len(nonzero_indices) > 0:
|
||||
plt.semilogy(nonzero_indices, np.abs(filtered_coeffs[nonzero_indices]), 'o-',
|
||||
label=f'Полином {degree}-й степени', markersize=6)
|
||||
|
||||
plt.axhline(y=MIN_COEFF_ABS, color='r', linestyle='--', alpha=0.7, label=f'Порог {MIN_COEFF_ABS:.0e}')
|
||||
plt.xlabel('Индекс коэффициента')
|
||||
plt.ylabel('Абсолютное значение коэффициента')
|
||||
plt.title('Коэффициенты после фильтрации (логарифмическая шкала)')
|
||||
plt.legend()
|
||||
plt.grid(True, alpha=0.3)
|
||||
|
||||
# График 5: Влияние фильтрации на ошибку
|
||||
plt.subplot(2, 3, 5)
|
||||
original_errors = []
|
||||
filtered_errors = []
|
||||
|
||||
for degree in degrees:
|
||||
poly_original = np.poly1d(coefficients[degree])
|
||||
poly_filtered = np.poly1d(filtered_coefficients[degree])
|
||||
|
||||
predicted_original = poly_original(adc_values)
|
||||
predicted_filtered = poly_filtered(adc_values)
|
||||
|
||||
error_original = np.max(np.abs(predicted_original - temperatures))
|
||||
error_filtered = np.max(np.abs(predicted_filtered - temperatures))
|
||||
|
||||
original_errors.append(error_original)
|
||||
filtered_errors.append(error_filtered)
|
||||
|
||||
x_pos = np.arange(len(degrees))
|
||||
width = 0.35
|
||||
|
||||
plt.bar(x_pos - width/2, original_errors, width, label='Оригинальные коэф.', alpha=0.7)
|
||||
plt.bar(x_pos + width/2, filtered_errors, width, label='После фильтрации', alpha=0.7)
|
||||
|
||||
plt.xlabel('Степень полинома')
|
||||
plt.ylabel('Максимальная ошибка (°C)')
|
||||
plt.title('Влияние фильтрации коэффициентов на точность')
|
||||
plt.xticks(x_pos, degrees)
|
||||
plt.legend()
|
||||
plt.grid(True, alpha=0.3)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
# Вывод численных результатов
|
||||
print("=" * 80)
|
||||
print("РЕЗУЛЬТАТЫ АППРОКСИМАЦИИ С ФИЛЬТРАЦИЕЙ КОЭФФИЦИЕНТОВ")
|
||||
print("=" * 80)
|
||||
|
||||
for degree in degrees:
|
||||
coeffs = coefficients[degree]
|
||||
filtered_coeffs = filtered_coefficients[degree]
|
||||
|
||||
poly_original = np.poly1d(coeffs)
|
||||
poly_filtered = np.poly1d(filtered_coeffs)
|
||||
|
||||
predicted_original = poly_original(adc_values)
|
||||
predicted_filtered = poly_filtered(adc_values)
|
||||
|
||||
max_error_original = np.max(np.abs(predicted_original - temperatures))
|
||||
rms_error_original = np.sqrt(np.mean((predicted_original - temperatures)**2))
|
||||
|
||||
max_error_filtered = np.max(np.abs(predicted_filtered - temperatures))
|
||||
rms_error_filtered = np.sqrt(np.mean((predicted_filtered - temperatures)**2))
|
||||
|
||||
# Подсчет нулевых коэффициентов
|
||||
zero_count = np.sum(filtered_coeffs == 0)
|
||||
total_count = len(filtered_coeffs)
|
||||
|
||||
print(f"\nПолином {degree}-й степени:")
|
||||
print(f"Максимальная ошибка: {max_error_original:.3f}°C -> {max_error_filtered:.3f}°C")
|
||||
print(f"Среднеквадратичная ошибка: {rms_error_original:.3f}°C -> {rms_error_filtered:.3f}°C")
|
||||
print(f"Обнулено коэффициентов: {zero_count}/{total_count}")
|
||||
|
||||
print("Коэффициенты после фильтрации:")
|
||||
for i, coeff in enumerate(filtered_coeffs):
|
||||
power = len(filtered_coeffs) - i - 1
|
||||
status = "✓" if coeff != 0 else "✗ (обнулен)"
|
||||
print(f" a_{power} = {coeff:.6e} {status}")
|
||||
|
||||
# Рекомендация
|
||||
print("\n" + "=" * 80)
|
||||
print("РЕКОМЕНДАЦИЯ:")
|
||||
print(f"Порог обнуления коэффициентов: {MIN_COEFF_ABS:.0e}")
|
||||
print("Полином 3-й степени обеспечивает оптимальный баланс между точностью")
|
||||
print("и сложностью реализации. Фильтрация малых коэффициентов практически")
|
||||
print("не влияет на точность, но упрощает реализацию на embedded системах.")
|
||||
print("=" * 80)
|
||||
11
Информация для программиста (УПП СП СЭД)/CALC/untitled.m
Normal file
11
Информация для программиста (УПП СП СЭД)/CALC/untitled.m
Normal file
@@ -0,0 +1,11 @@
|
||||
clc
|
||||
clear all
|
||||
|
||||
|
||||
alpha = 0.01
|
||||
ts = 500;
|
||||
|
||||
tau = (ts/1000000) * (1-alpha)/alpha
|
||||
|
||||
|
||||
alpha2 = ts/1000000 / (ts/1000000 + tau )
|
||||
Reference in New Issue
Block a user