Фильтр Баттерворта

@Arhel · Регистрация: 04.01.2012

Студворк — интернет-сервис помощи студентам

Доброго времени суток!

Подскажите, пожалуйста, как рассчитать коэффициенты полосового фильтра Баттерворта 8-го порядка методом билинейной трансформации? (Частота дискретизации = 3000 Гц, f1 = 450Гц, f2 = 750 Гц)

Пробовал использовать функцию fdatool, но, к сожалению, там нельзя выбрать какой-либо метод - билинейной трансформации или к примеру импульсной инвариации.

Можно, конечно, вручную..но лень =)

Спасибо!

@Arhel · 08.01.2012, 13:17 **[ТС]**

Друзья, подскажите, пожалуйста как восстановить аудиосигнал с помощью, например, фильтра Баттверворта? Составил код, который обрезает высокочастотную составляющую, но почему-то шум все равно присутствует в сигнале =(

Matlab M
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
[x,fs,nbits] = wavread('C:\sample3a_sk1352.wav');
x = x(:,1); % Select left channel only
nsamples = length(x);
 
xss = x(1:5:nsamples); % Sub-sample
fss = fs / 5; % Corresponding sampling frequency
 
f=21000; % 21KHz
x = 0.9*x + 0.1*cos(2*pi*f*(1:nsamples)'/fs);
 
wavplay(x, fs);
 
xss = x(1:5:nsamples);
fss = fs / 5; % Corresponding sampling frequency
wavplay(xss, fss);
 
[b,a] = butter(10,1/5);
xf = filter(b,a,x);
xfss = xf(1:5:nsamples);
fss = fs / 5; % Corresponding sampling frequency
wavplay(xfss, fss);
 
y = zeros(nsamples,1);
y(1:5:nsamples) = xfss;
wavplay(y, fs);
 
yf = filter(5*b,a,y);
wavplay(yf, fs);

В приложении зв. файл

@vital792 · 08.01.2012, 13:39

Подскажите, пожалуйста, как рассчитать коэффициенты полосового фильтра Баттерворта 8-го порядка методом билинейной трансформации? (Частота дискретизации = 3000 Гц, f1 = 450Гц, f2 = 750 Гц)

я сделал бы как то так:

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Fs = 3000;
t = 0 : 1/Fs : 4095/Fs;
signal = sin(2*pi*100*t)+sin(2*pi*300*t)+sin(2*pi*600*t)+0.5*randn(1, length(t));
% это типа сигнал=)
spec = fft(signal);
plot(abs(spec(1:floor(length(spec)/2))));
 
n = 8;
Wn = [450 750]/Fs;
[b, a] = butter(n, Wn, 'bandpass');
 
filtredSignal = filter(b, a, signal);
 
specFiltred = fft(filtredSignal);
hold on;
plot(abs(specFiltred(1:floor(length(specFiltred)/2))), 'r');

функция butter тоже не позволяет выбирать метод синтеза фильтра и использует как раз билинейное z преобразование. Если надо другим методом - придется сначала создать аналоговый прототип, и дискретизировать его тем же инвариантным преобразованием ИХ.

Можно, конечно, вручную..но лень =)

да и зачем? велосипед уже изобрели)

@Arhel · 08.01.2012, 13:47 **[ТС]**

Ой, спасибо огромное =)
Вообще да, мне еще надо методом импульсной инвариации...но сначала разобраться б с билинейным)
а как из этого получить коэффициенты фильтра?

Да уж..хотя я тут потренировался и в ручную рассчитал FIR bandpass 32 порядка =)хватило

@vital792 · 08.01.2012, 14:05

Сообщение от Arhel

xss = x(1:5:nsamples); % Sub-sample
fss = fs / 5; % Corresponding sampling frequency

а для чего все эти трюки?

Сообщение от Arhel

f=21000; % 21KHz
x = 0.9*x + 0.1*cos(2*pi*f*(1:nsamples)'/fs);

то есть вы добавляете к сигналу тональную помеху частотой 21 кГц

Сообщение от Arhel

[b,a] = butter(10,1/5);
xf = filter(b,a,x);

затем фильтруете сигнал. Все. ваша помеха исчезла - сравните спектр x и xf. А на слух, вы слышите частоту 21 кГц? Я лично нет. Мой порог около 17-18, дольше я не слышу

@vital792 · 08.01.2012, 14:11

попробовал, фильтр получился неустойчивым. За это я ненавижу рекурсивные фильтры.

Сообщение от Arhel

ручную рассчитал FIR bandpass 32 порядка =)

если можно использовать fir - это надо делать!!) Я вообще в своей практике никогда не сталкивался с задачами, когда fir фильтр не справлялся

@Arhel · 08.01.2012, 14:14 **[ТС]**

Да, сравнил спектры, помехи нет..но при прослушивании wav файла все равно звучит высокочастотный шум =( а как вывести коэффициенты фильтра?

@Arhel · 08.01.2012, 14:15 **[ТС]**

=)буду знать..только начинаю пользоваться матлабом..вот и столько вопросов) не подскажешь, как вывести коэффициенты фильтра?

@vital792 · 08.01.2012, 14:23

Сообщение от Arhel

а как вывести коэффициенты фильтра?

так ведь [b, a] это они и есть. Числитель и знаменатель.

@Arhel · 08.01.2012, 19:31 **[ТС]**

Да,а как их значения вывести на экран? =)

@vital792 · 08.01.2012, 19:39

Сообщение от Arhel

Да,а как их значения вывести на экран? =)

вот так:

Code
1
b, a

)) Выведутся в command window
Или в каком виде вывести?

@Arhel · 08.01.2012, 20:33 **[ТС]**

а,да..что-то я уже запутался..
спасибо большое ;-)

Добавлено через 45 минут
Извини,что столько много вопросов..
У меня получились коэффициенты:
b =

1.0e-004 *

Columns 1 through 10

0.0018 0 -0.0141 0 0.0494 0 -0.0987 0 0.1234 0

Columns 11 through 17

-0.0987 0 0.0494 0 -0.0141 0 0.0018

a =

1.0e+003 *

Columns 1 through 10

0.0010 -0.0118 0.0673 -0.2460 0.6444 -1.2801 1.9930 -2.4792 2.4896 -2.0247

Columns 11 through 17

1.3292 -0.6973 0.2867 -0.0894 0.0200 -0.0029 0.0002

А порядок фильтра = 8. По идее b, a должны быть по 9? получается что каждого по 17...

@vital792 · 08.01.2012, 22:04

Сообщение от Arhel

А порядок фильтра = 8. По идее b, a должны быть по 9? получается что каждого по 17...

Верное замечание - а я и не обратил внимание. Полосовой фильтр расчитывается порядка 2*N

% If Wn is a two-element vector, Wn = [W1 W2], BUTTER returns an
% order 2N bandpass filter with passband W1 < W < W2.

@Arhel · 08.01.2012, 22:27 **[ТС]**

Да,все получилось! =)ееее))пасиб

@Arhel · 11.01.2012, 21:57 **[ТС]**

Vital792,
снова привет!подскажи,плиз, правильно ли я подправил метод билинейного преобразования на метод импульсной инвариантности?и еще вопрос..сигнал ты откуда брал?

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Fs = 3000;
t = 0 : 1/Fs : 4095/Fs;
signal = sin(2*pi*100*t)+sin(2*pi*300*t)+sin(2*pi*600*t)+0.5*randn(1, length(t));[/code]
% это типа сигнал=)
[code]spec = fft(signal);
plot(abs(spec(1:floor(length(spec)/2))));
 
n = 4;
Wn = [450 750]/Fs;
[b, a] = butter(n, Wn, 'bandpass');
[bz,az] = impinvar(b,a,Fs); 
 
filtredSignal = filter(b, a, signal);
 
specFiltred = fft(filtredSignal);
hold on;
plot(abs(specFiltred(1:floor(length(specFiltred)/2))), 'r');

@vital792 · 11.01.2012, 22:08

Сообщение от Arhel

подскажи,плиз, правильно ли я подправил метод билинейного преобразования на метод импульсной инвариантности?

как сказал Арнольд Шварцнегер в фильме "Тегминатор":
"НЕТ"

Сообщение от Arhel

и еще вопрос..сигнал ты откуда брал?

так ведь и не от куда) я же его генерирую:

Сообщение от Arhel

t = 0 : 1/Fs : 4095/Fs;

задаю время

Сообщение от Arhel

signal = sin(2*pi*100*t)+sin(2*pi*300*t)+sin(2*pi *600*t)+0.5*randn(1, length(t));

генерю сумму из трех синусоид с разными частотами + белый шум. Это и есть сигнал. На нем и проверял

@Arhel · 11.01.2012, 22:12 **[ТС]**

Понятно..пасиб)а что не так в коде?вроде бы использую команду impinvar

@vital792 · 11.01.2012, 22:29

Сообщение от Arhel

вроде бы использую команду impinvar

команда дискретизирует аналоговый фильтр, а ты получаешь уже готовые коэффициенты цифрового. Зайди в код функции butter, и посмотри как она работает. Обрати внимание на параметр 's' Извини, щас под рукой нет матлаба, ниче по проверить не могу( Тока завтра

@Arhel · 11.01.2012, 22:45 **[ТС]**

Ок, спасибо..буду вникать)

@Arhel · 19.01.2012, 01:47 **[ТС]**

Снова привет!Не мог бы проверить следующий код? Методом инвариантности получается все-таки какая-то фигня..не могу понять где ошибка =( фильтр получается совершенно иным, чем при билинейной трансформации..

Matlab M
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Fs = 3000;
t = 0 : 1/Fs : 4095/Fs;
signal = sin(2*pi*200*t)+sin(2*pi*500*t)+sin(2*pi*900*t)+0.5*randn(1, length(t));
% это типа сигнал=)
spec = fft(signal);
plot(abs(spec(1:floor(length(spec)/2))));
 
n = 4;
Wn = [450 750]/Fs;
[b,a] = butter(n, Wn, 'bandpass','s');
[bz,az] = impinvar(b,a,Fs); 
 
 
filtredSignal = filter(bz, az, signal);
 
specFiltred = fft(filtredSignal);
hold on;
plot(abs(specFiltred(1:floor(length(specFiltred)/2))), 'r');
fvtool(b,a)

Новые блоги и статьи Все статьи Все блоги /
SDL3 для Web (WebAssembly): Загрузка PNG с прозрачным фоном с помощью SDL3_image 8Observer8 10.02.2026 Содержание блога Библиотека SDL3_image содержит инструменты для расширенной работы с изображениями. Пошагово создадим проект для загрузки изображения формата PNG с альфа-каналом (с прозрачным. . .	Установка Qt-версии Lazarus IDE в Debian Trixie Xfce volvo 10.02.2026 В общем, достали меня глюки IDE Лазаруса, собранной с использованием набора виджетов Gtk2 (конкретно: если набирать текст в редакторе и вызвать подсказку через Ctrl+Space, то после закрытия окошка. . .	SDL3 для Web (WebAssembly): Работа со звуком через SDL3_mixer 8Observer8 08.02.2026 Содержание блога Пошагово создадим проект для загрузки звукового файла и воспроизведения звука с помощью библиотеки SDL3_mixer. Звук будет воспроизводиться по клику мышки по холсту на Desktop и по. . .	SDL3 для Web (WebAssembly): Основы отладки веб-приложений на SDL3 по USB и Wi-Fi, запущенных в браузере мобильных устройств 8Observer8 07.02.2026 Содержание блога Браузер Chrome имеет средства для отладки мобильных веб-приложений по USB. В этой пошаговой инструкции ограничимся работой с консолью. Вывод в консоль - это часть процесса. . .
SDL3 для Web (WebAssembly): Обработчик клика мыши в браузере ПК и касания экрана в браузере на мобильном устройстве 8Observer8 02.02.2026 Содержание блога Для начала пошагово создадим рабочий пример для подготовки к экспериментам в браузере ПК и в браузере мобильного устройства. Потом напишем обработчик клика мыши и обработчик. . .	Философия технологии iceja 01.02.2026 На мой взгляд у человека в технических проектах остается роль генерального директора. Все остальное нейронки делают уже лучше человека. Они не могут нести предпринимательские риски, не могут. . .	SDL3 для Web (WebAssembly): Вывод текста со шрифтом TTF с помощью SDL3_ttf 8Observer8 01.02.2026 Содержание блога В этой пошаговой инструкции создадим с нуля веб-приложение, которое выводит текст в окне браузера. Запустим на Android на локальном сервере. Загрузим Release на бесплатный. . .	SDL3 для Web (WebAssembly): Сборка C/C++ проекта из консоли 8Observer8 30.01.2026 Содержание блога Если вы откроете примеры для начинающих на официальном репозитории SDL3 в папке: examples, то вы увидите, что все примеры используют следующие четыре обязательные функции, а. . .

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15

	Фильтр Баттерворта 07.01.2012, 23:17. Показов 19836. Ответов 21 Метки нет (Все метки) Доброго времени суток! Подскажите, пожалуйста, как рассчитать коэффициенты полосового фильтра Баттерворта 8-го порядка методом билинейной трансформации? (Частота дискретизации = 3000 Гц, f1 = 450Гц, f2 = 750 Гц) Пробовал использовать функцию fdatool, но, к сожалению, там нельзя выбрать какой-либо метод - билинейной трансформации или к примеру импульсной инвариации. Можно, конечно, вручную..но лень =) Спасибо! 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	08.01.2012, 13:47 [ТС]
	Ой, спасибо огромное =) Вообще да, мне еще надо методом импульсной инвариации...но сначала разобраться б с билинейным) а как из этого получить коэффициенты фильтра? Да уж..хотя я тут потренировался и в ручную рассчитал FIR bandpass 32 порядка =)хватило 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	08.01.2012, 14:14 [ТС]
	Да, сравнил спектры, помехи нет..но при прослушивании wav файла все равно звучит высокочастотный шум =( а как вывести коэффициенты фильтра? 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	08.01.2012, 14:15 [ТС]
	=)буду знать..только начинаю пользоваться матлабом..вот и столько вопросов) не подскажешь, как вывести коэффициенты фильтра? 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	08.01.2012, 19:31 [ТС]
	Да,а как их значения вывести на экран? =) 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	08.01.2012, 20:33 [ТС]
	а,да..что-то я уже запутался.. спасибо большое ;-) Добавлено через 45 минут Извини,что столько много вопросов.. У меня получились коэффициенты: b = 1.0e-004 * Columns 1 through 10 0.0018 0 -0.0141 0 0.0494 0 -0.0987 0 0.1234 0 Columns 11 through 17 -0.0987 0 0.0494 0 -0.0141 0 0.0018 a = 1.0e+003 * Columns 1 through 10 0.0010 -0.0118 0.0673 -0.2460 0.6444 -1.2801 1.9930 -2.4792 2.4896 -2.0247 Columns 11 through 17 1.3292 -0.6973 0.2867 -0.0894 0.0200 -0.0029 0.0002 А порядок фильтра = 8. По идее b, a должны быть по 9? получается что каждого по 17... 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	08.01.2012, 22:27 [ТС]
	Да,все получилось! =)ееее))пасиб 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	11.01.2012, 22:12 [ТС]
	Понятно..пасиб)а что не так в коде?вроде бы использую команду impinvar 0

@Arhel 0 / 0 / 0 Регистрация: 04.01.2012 Сообщений: 15
	11.01.2012, 22:45 [ТС]
	Ок, спасибо..буду вникать) 0