Как вручную построить частотную характеристику полосового фильтра Баттерворта в MATLAB без функции freqz?

15

У меня есть код, как показано ниже, который применяет полосовой фильтр к сигналу. Я довольно новичок в DSP, и я хочу понять, что происходит за кулисами, прежде чем я продолжу.

Для этого я хочу знать, как построить частотную характеристику фильтра без использования freqz.

[b, a] = butter(order, [flo fhi]);
filtered_signal = filter(b, a, unfiltered_signal)

Учитывая результаты, [b, a]как бы я это сделал? Кажется, это будет простая задача, но мне трудно найти то, что мне нужно, в документации или онлайн.

Я также хотел бы понять, как сделать это как можно быстрее, например, с помощью fftили другого быстрого алгоритма.

Уильям
источник

Ответы:

25

Мы знаем, что в целом передаточная функция фильтра определяется как:

H(z)=k=0Mbkzkk=0Nakzk

Теперь подставим z=ejω чтобы оценить передаточную функцию на единичном круге:

ЧАС(еJω)знак равноΣКзнак равно0MбКе-JωКΣКзнак равно0NaКе-JωК

Таким образом, это становится только проблемой полиномиальной оценки при заданном ω . Вот шаги:

  • Создайте вектор угловых частот ωзнак равно[0,...,π] для первой половины спектра (не нужно подниматься до 2π ) и сохраните его в w.
  • Предварительно вычислите показатель степени е-Jω для всех из них и сохраните его в переменной ze.
  • Используйте polyvalфункцию для вычисления значений числителя и знаменателя, вызывая:, polyval(b, ze)разделите их и сохраните в H. Поскольку нас интересует амплитуда, то берем абсолютное значение результата.
  • Преобразовать в шкалу дБ с помощью: ЧАСdВзнак равно20журнал10ЧАС - в этом случае 1 является эталонным значением.

Поместить все это в код:

%% Filter definition
a = [1 -0.5 -0.25]; % Some filter with lot's of static gain
b = [1 3 2];

%% My freqz calculation
N = 1024; % Number of points to evaluate at
upp = pi; % Evaluate only up to fs/2
% Create the vector of angular frequencies at one more point.
% After that remove the last element (Nyquist frequency)
w = linspace(0, pi, N+1); 
w(end) = [];
ze = exp(-1j*w); % Pre-compute exponent
H = polyval(b, ze)./polyval(a, ze); % Evaluate transfer function and take the amplitude
Ha = abs(H);
Hdb  = 20*log10(Ha); % Convert to dB scale
wn   = w/pi;
% Plot and set axis limits
xlim = ([0 1]);
plot(wn, Hdb)
grid on

%% MATLAB freqz
figure
freqz(b,a)

Исходный вывод freqz:

введите описание изображения здесь

И вывод моего скрипта:

введите описание изображения здесь

И быстрое сравнение в линейном масштабе - выглядит великолепно!

[h_f, w_f] = freqz(b,a);
figure
xlim = ([0 1]);
plot(w, Ha) % mine
grid on
hold on
plot(w_f, abs(h_f), '--r') % MATLAB
legend({'my freqz','MATLAB freqz'})

введите описание изображения здесь

Теперь вы можете переписать его в некоторую функцию и добавить несколько условий, чтобы сделать его более полезным.


Другой способ (ранее предложенный является более надежным) заключается в использовании фундаментального свойства, заключающегося в том, что частотная характеристика фильтра является преобразованием Фурье его импульсной характеристики:

H(ω)=F{h(t)}

Поэтому вы должны подать в вашу систему сигнал δ(t) , рассчитать отклик вашего фильтра и взять его БПФ:

d = [zeros(1,length(w_f)) 1 zeros(1,length(w_f)-1)];
h = filter(b, a, d);
HH = abs(fft(h));
HH = HH(1:length(w_f));

Для сравнения это даст следующее:

введите описание изображения здесь

Jojek
источник
1
Подробное объяснение
Партида
Я думаю об этой линии a = [1 -0.5 -0.25]; % Some filter with lot's of static gain. Можете ли вы объяснить мне выбор этих параметров здесь, пожалуйста. Я читаю руководство по моему Matlab, и оно написано [h,w] = freqz(hfilt,n)в одной части синапса. Вы даете два фильтра (а, б) в freqz. Оба фильтра в hfilt? Или один в n?
Лео Леопольд Герц 준영
Просто чтобы уточнить для других: «Не нужно идти до 2 пи», когда коэффициенты реальны. Существуют приложения для фильтров со сложными коэффициентами, и в этом случае спектр больше не будет симметричным, и в этом случае потребуется увеличить частоту до 2 пи.
Дэн Бошен
14

это всего лишь дополнение к ответу Джохека, которое является более общим и совершенно хорошим, когда используется математика двойной точности. когда точность меньше, возникает «проблема косинуса», которая возникает, когда либо частота в частотной характеристике очень низкая (намного ниже, чем Найквист), а также когда резонансные частоты фильтра очень низкие.

|ЧАС(еJω)|2|ЧАС(е-Jω)|знак равно|ЧАС(еJω)|

рассмотрим этот триг личности:

соз(ω) знак равно 1-2грех2(ω2)

грех2(ω2)ω0

грех2(ω2)

ЧАС(Z)знак равноб0+б1Z-1+б2Z-2a0+a1Z-1+a2Z-2

который имеет сложную частотную характеристику

ЧАС(еJω)знак равноб0+б1е-Jω+б2е-J2ωa0+a1е-Jω+a2е-J2ω

который имеет величину в квадрате:

|ЧАС(еJω)|2знак равно|б0+б1е-Jω+б2е-J2ω|2|a0+a1е-Jω+a2е-J2ω|2знак равно(б0+б1соз(ω)+б2соз(2ω))2+(б1грех(ω)+б2грех(2ω))2(a0+a1соз(ω)+a2соз(2ω))2+(a1грех(ω)+a2грех(2ω))2знак равноб02+б12+б22+2б1(б0+б2)соз(ω)+2б0б2соз(2ω)a02+a12+a22+2a1(a0+a2)соз(ω)+2a0a2соз(2ω)

|ЧАС(еJω)|соз(ω)соз(2ω)ω11

используя приведенную выше идентичность триггера, вы получаете для квадрата величины:

|H(ejω)|2=b02+b12+b22+2b1(b0+b2)cos(ω)+2b0b2cos(2ω)a02+a12+a22+2a1(a0+a2)cos(ω)+2a0a2cos(2ω)=b02+b12+b22+2b1(b0+b2)(12sin2(ω2))+2b0b2(12sin2(ω))a02+a12+a22+2a1(a0+a2)(12sin2(ω2))+2a0a2(12sin2(ω))=b02+b12+b22+2b1(b0+b2)(12sin2(ω2))+2b0b2(2cos2(ω)1)a02+a12+a22+2a1(a0+a2)(12sin2(ω2))+2a0a2(2cos2(ω)1)=b02+b12+b22+2b1(b0+b2)(12sin2(ω2))+2b0b2(2(12sin2(ω2))21)a02+a12+a22+2a1(a0+a2)(12sin2(ω2))+2a0a2(2(12sin2(ω2))21)=b02+b12+b22+2b1(b0+b2)(12ϕ)+2b0b2(2(12ϕ)21)a02+a12+a22+2a1(a0+a2)(12ϕ)+2a0a2(2(12ϕ)21)=b02+b12+b22+2b1(b0+b2)(12ϕ)+2b0b2(18ϕ+8ϕ2)a02+a12+a22+2a1(a0+a2)(12ϕ)+2a0a2(18ϕ+8ϕ2)=b02+b12+b22+2b1b0+2b1b24(b1b0+b1b2)ϕ+2b0b216b0b2ϕ+16b0b2ϕ2a02+a12+a22+2a1a0+2a1a24(a1a0+a1a2)ϕ+2a0a216a0a2ϕ+16a0a2ϕ2=(b02+b12+b22+2b1b0+2b1b2+2b0b2)4(b1b0+b1b24b0b2)ϕ+16b0b2ϕ2(a02+a12+a22+2a1a0+2a1a2+2a0a2)4(a1a0+a1a24a0a2)ϕ+16a0a2ϕ2=14(b02+b12+b22+2b1b0+2b1b2+2b0b2)(b1b0+b1b24b0b2)ϕ+4b0b2ϕ214(a02+a12+a22+2a1a0+2a1a2+2a0a2)(a1a0+a1a24a0a2)ϕ+4a0a2ϕ2=(b0+b1+b22)2ϕ(4b0b2(1ϕ)+b1(b0+b2))(a0+a1+a22)2ϕ(4a0a2(1ϕ)+a1(a0+a2))

where ϕsin2(ω2)

if your gear is intending to plot this as dB, it comes out as

20log10|H(ejω)| = 10log10((b0+b1+b22)2ϕ(4b0b2(1ϕ)+b1(b0+b2)))10log10((a0+a1+a22)2ϕ(4a0a2(1ϕ)+a1(a0+a2)))

so your division turns into subtraction, but you have to be able to compute logarithms to some base or another. numerically, you will have much less trouble with this for low frequencies than doing it the apparent way.

robert bristow-johnson
источник
2
That's really cool, thank you Robert! +1
jojek
@Robert I "believe" similar to my comment for Jojek above that this only applies as well when the coefficients are real (and therefore the spectrum is symmetric and thus the magnitude converts to cosines as you show)... Am I correct?
Dan Boschen
yes. that commitment is made when you go from the first line of |H(ejω)|2=... to the second line. no going back after that.
Роберт Бристоу-Джонсон