UPP/Информация для программиста (УПП СП СЭД)/CALC/BiquadFilterDesigner.m

%% КЛАСС ДЛЯ РАСЧЕТА КОЭФФИЦИЕНТОВ БИКВАДРАТНЫХ ФИЛЬТРОВ
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