Цифровой фильтр Баттерворта

[Devil]

Здравствуйте! Синтезирую цифровой фильтр. В процессе расчета фильтра Баттерворта 2 порядка у меня возникли следующие вопросы:
1. меняются ли коэффициенты передаточной функции в пределах одного порядка?
2. возможен ли переход от аналогового фильтра к цифровому, минуя билинейное преобразование?
3. чем может быть обусловлен подъем ачх в полосе пропускания, если фильтр обычный (без усиления)?
Буду благодарна за любую информацию.

[Бахурин Сергей]

Здравствуйте! Синтезирую цифровой фильтр. В процессе расчета фильтра Баттерворта 2 порядка у меня возникли следующие вопросы:
1. меняются ли коэффициенты передаточной функции в пределах одного порядка?

не совсем понял что имеется ввиду? Если заданы одни и те же исходные данные то коэффициенты фильтра не меняются. Если меняете что-то в исходных данных то и к-ты поменяются тоже.

2. возможен ли переход от аналогового фильтра к цифровому, минуя билинейное преобразование?

Есть еще метод расчета на основе инвариантной импульсной характеристики, но он редко применяется сегодня поскольку билинейное преобразование отлично реализуется программно.

3. чем может быть обусловлен подъем ачх в полосе пропускания, если фильтр обычный (без усиления)?
Буду благодарна за любую информацию.

Многими факторами: неправильный расчет к-тов фильтра, неправильная нормировка к-тов фильтра и т.д. надо смотреть в каждом отдельном случае.

[Devil]

Спасибо Вам, Сергей. Все получилось, очень помогли статьи на сайте.

[pvr]

Здравствуйте!
Мне нужно реализовать цифровой фильтр Баттерворта 6-го порядка с частотой среза fp на уровне 3 дБ, равной 35 Гц. Частота дискретизации fs = 256 Гц.
Не получается.
Давайте, чтобы было проще, реализуем то же самое, но для фильтра 3-порядка.
По условию epsp = 1, откуда alpha тоже равна 1. Получаем
Hн(s) = 1 / (1 + 2 * s + 2 * s^2 + s^3).
Теперь приведём частоту среза:
Fp = fp / (fs / 2) = 0.2734375.
Сделаем подстановку s -> s / Fp, получим
H = 1 / (1 + 7.3142 * s + 26.7493 * s^2 + 48.9131 * s^3).
Теперь, если сделать билинейное преобразование
s = (2 / T) * (1 - z) / (1 + z),
то получится
H(z) = (1 + 3 * z + 3 * z^2 + 1 * z^3) /
(6572029952 - 19702050816 * s + 19688026112 * s^2 - 6558005760 * s^3).
После приведения к а0 = 1 в числителе получаются практически нули.

Что я делаю не так?

[pvr]

Фильтр правильно срезает, если в билинейном преобразовании выбрать T ~ 3 / 8. Ума не приложу, почему это так...

[Бахурин Сергей]

вот матлаб пример который рассчитывает ваш фильтр в точности как в примере. Пошагово сможете отладить ваш расчет.

Код: Выделить всё

clear all; close all; clc;
%% Исходные данные
% Данный пример расчитывает ФНЧ согласно примеру
% http://www.dsplib.ru/content/filters/butterex/butterex.html
% Для других исходных данных
wp = 35;    % Частота среза (Гц)
Fs = 256;   % Частота дискретизации
Rp = 3;     % Неравномерность АЧХ дБ
N  = 5;     % Порядок фильтра


%% Расчет нормированного анлогового ФНЧ
r = mod(N,2);
L = floor(N/2);

ep = sqrt(10^(Rp/10) - 1);
alpha = 1/(ep^(1/N));

n = 1:L;
teta = pi*(2*n-1)/(2*N);

if r==1
    a = ep*[1 alpha];
else
    a = ep;
end
for n = 1:L
   a = conv(a, [1 2*alpha*sin(teta(n)) alpha^2]); 
end
a = fliplr(a);
b = zeros(1,length(a));
b(1) = 1;


%% Пересчет частоты и преобразование ФНЧ-ФНЧ
wp = 2*wp/Fs;
Wp = tan(pi*wp/2);
[b, a] = subst(b, a, [0 1], [Wp 0])

%% Билинейное преобразование
[b, a] = subst(b, a, [1 -1], [1 1]);

%% построение характеристики фильтра
fvtool(b,a,'Fs', Fs)

функция рациональной подстановки (смотри подробнее здесь)

Код: Выделить всё

function [bs, as] = subst(b, a, beta, alpha)
% Rational substitution.
% [BS, AS] = SUBST(B, A, BETA, ALPHA) 
% function calculates the following substitution
%
%                       N
%                      SUM  b(n) * R(z)^n
%                      n=0
%   G(z) = H(R(z)) = ----------------------- , here
%                       N
%                      SUM  a(m) * R(z)^m
%                      m=0
%
%             N                                  P
%            SUM  b(n) * s^n                    SUM  beta(p) * z^p
%            n=0                                p=0
%   H(s) = ----------------------- ,     R(z) = -------------------
%             M                                  Q
%            SUM  a(m) * s^m                    SUM alpha(q) * z^q
%            m=0                                q=0
%
% B     - coefficients vector of H(s) numerator  
% A     - coefficients vector of H(s) denominator 
% BETA  - coefficients vector of R(z) numerator 
% ALPHA - coefficients vector of R(z) denominator
%
% BS    -  coefficients vector of H(R(z)) numerator  (the result)
% AS    -  coefficients vector of H(R(z)) denominator (the result)
% ------------------------------------------------------------------------
%
% Дробно-рациональная подстановка
% [BS, AS] = SUBST(B, A, BETA, ALPHA) расчитывает подстановку вида
%
%                       N
%                      SUM  b(n) * R(z)^n
%                      n=0
%   G(z) = H(R(z)) = ----------------------- , где
%                       N
%                      SUM  a(m) * R(z)^m
%                      m=0
%
%             N                                  P
%            SUM  b(n) * s^n                    SUM  beta(p) * z^p
%            n=0                                p=0
%   H(s) = ----------------------- , a   R(z) = -------------------
%             M                                  Q
%            SUM  a(m) * s^m                    SUM alpha(q) * z^q
%            m=0                                q=0
% ------------------------------------------------------------------------
% www.dsplib.ru
%


LA = length(alpha);
LB = length(beta);

if(LA > LB)
    beta = [beta, zeros(1, LA - LB)];
elseif(LB > LA)
    alpha = [alpha zeros(1, LB - LA)];
    LA = LB;
end

N = length(a) ;
M = N * LA - N - LA + 2;

NUM = zeros(N, M);
NUM(1, 1) = 1;
DEN = NUM;

for k =  2 : N
   x = conv(NUM(k - 1, :), beta);
   NUM(k, :) = x(1 : M);
   x = conv(DEN(k - 1, :), alpha);
   DEN(k,:) = x(1 : M);
end

ND = zeros(N, M);
for k = 1 : N
   x = conv(DEN(N - k + 1, :), NUM(k, :)); 
   ND(k, :) = x(1 : M);    
end

as = a * ND;
bs = b * ND;

Получился такой фильтр: