В модель добавлена библиотека CMSIS-DSP и вообще все либы CMSIS

This commit is contained in:
2025-11-13 17:14:43 +03:00
parent 75bed20511
commit 5299cc5b12
1074 changed files with 380657 additions and 7 deletions

View File

@@ -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

View 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);