Разработка КИХ–фильтра с помощью MATLAB

Метод взвешивания

Этапы вычисления коэффициентов стандартного частотно-избирательного КИХ–фильтра с линейной фазовой характеристикой с помощью метода взвешивания можно упорядочить следующим образом.

1) Задать необходимую частотную характеристику.

2) Выбрать весовую функцию и оценить число коэффициентов фильтра, N.

3) Получить идеальную частотную характеристику, h_D (n) (усеченную до N значений).

4) Получить N коэффициентов весовой функции, w(n).

5) Получить коэффициенты КИХ–фильтра, воздействовав на частотную характеристику весовой функцией, h(n) = h_D (n) * w(n).

Для стандартного частотно-избирательного КИХ–фильтра с линейной фазовой характеристикой (нижних частот, верхних частот, полосовые и режекторные фильтры), при разработке которого используется метод взвешивания, ключевой командой высокого уровня в Toolbox является команда fir1. Синтаксис данной команды (в стандартной форме):

b = fir1(N-l,Fc).

Команда в стандартной форме вычисляет и возвращает коэффициенты N-точечной импульсной характеристики КИХ–фильтра с частотой среза Fc. Команда возвращает коэффициенты в вектор b, упорядоченный по возрастанию отрицательной степени z:

Параметр команды N–1 задает порядок фильтра (обычно на единицу меньше числа коэффициентов КИХ–фильтра). Частота среза Fc нормирована на частоту Найквиста (т.е. половину частоты дискретизации).

По умолчанию стандартная команда fir1 действует на данные весовой функцией Хэмминга и в ней предполагается использование фильтра нижних частот (или полосовой фильтр, если Fc задает более одной частоты среза). Стандартную команду можно расширить, задав тип фильтра и/или весовую функцию. В этом случае используется такой синтаксис:

b = fir1(N - 1, Fc, 'тип фильтра')

b = fir1(N - 1, Fc, window)

b = fir1(N - 1, Fc, 'тип фильтра', window)

Для фильтра верхних частот (highpass) тип задается словом 'high', для режекторного (stopband) – словом 'stop'. Для полосовых и режекторных фильтров переменная Fc – вектор, который задает частоты среза. Для фильтров верхних частот и режекторных длина фильтра должна быть нечетным целым числом (четные целые не подходят, поскольку, как показано ниже, это приведет к нулевой амплитудной характеристике на частоте Найквиста).

MATLAB поддерживает использование различных весовых функций, включая функции Хэмминга (Hamming), Хеннинга (Hanning), прямоугольную (boxcar), Кайзера (Kaiser) и Чебышева (Chebyshev). Для получения весовых коэффициентов следует использовать такой синтаксис:

w = boxcar(N)

w = blackman(N)

w = hamming(N)

w = hanning(N)

w = kaiser(N, beta)

На практике команда взвешивания часто вкладывается в команду fir1(см. примеры ниже).

Нужно отметить, что из-за различий в реализациях результаты, полученные при разработке КИХ–фильтров на основе методов вырезания с помощью MATLAB, могут отличаться от результатов, полученных с помощью других программ. Например, в MATLAB после вырезания коэффициенты импульсной характеристики могут масштабироваться для получения в середине полосы пропускания единичной амплитудно–частотной характеристики. Чтобы запретить такое поведение, следует добавить слово 'noscale', например:

b = firl(N - 1, Fc, 'noscale').

Кроме того, реализация в MATLAB большинства весовых функций может немного меняться, приводя к различным результатам. Разработчик должен знать об этих отличиях и при необходимости делать соответствующие допущения, чтобы их компенсировать.

Пример 1.Определите коэффициенты КИХ–фильтра нижних частот с линейной фазовой характеристикой с краями полосы пропускания (f_s) и полосы подавления (f_p) на частотах 1 и 4,3 кГц соответственно. Используйте весовую функцию Хэмминга, частоту дискретизации предполагайте равной 10 кГц.

Решение

Из табл. 3 находим подходящее соотношение между шириной перехода и длиной фильтра N для фильтра на основе весовой функции Хэмминга:

В нашем случае равно 0.33 (из (f_p- f_s)/F_s), так что длина фильтра N=10. Используя тот же подход, что и в теоретическом материале, принимаем, что реальная частота среза (с учетом эффекта размывания) лежит посередине между заданными краями полос пропускания и подавления, т.е. на частоте 2,65 кГц. В MATLAB частоту среза можно нормировать на половину частоты дискретизации. Тогда:

fc(нормированная)=2,65/5=0,53.

Команды MATLAB:

fc=0.53; % Частота среза (нормированная на Fs/2)

N= 10; % Длина фильтра

hd=fir1(N-1,fc,boxcar(N)); % Усеченная идеальная импульсная характеристика

wn=hamming(N); % Вычислить коэффициенты функции Хэмминга

hn=fir1(N-1,fc,wn); % Получить коэффициенты взвешенной функции

Пример 2. Определите коэффициенты и изобразите амплитудно-частотную характеристику полосового КИХ-фильтра, используя весовую функцию Кайзера. Фильтр должен удовлетворять спецификациям:

полоса пропускания 150-200 Гц

ширина полосы перехода 50 Гц

неравномерность в полосе пропускания 0.1 дБ

затухание в полосе подавления 60 дБ

частота дискретизации 1 кГц

Решение

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

дБ, откуда = 0.0115

дБ, откуда = 0.001

Следовательно,

Требования к затуханию удовлетворим при помощи функции Кайзера. Число коэффициентов фильтра в данном случае равно:

Пусть N=73. Параметр неравномерности выражается следующим образом:

Значения N и используются для работы функции kaiser, на выходе которой получаются значения w(n). Чтобы учесть эффект смазывания, при вычислении идеальной импульсной характеристики использовались частоты среза и , т.е =125 Гц и =275 Гц соответственно.

Следует отметить, что в данном примере коэффициенты весовой функции не вычисляются отдельно, просто в команду fir1 включен тип весовой функции.

Команды MATLAB:

FS=1000; %Частота дискретизации

FN=FS/2; %Частота Найквиста

N=73; %Длина фильтра

beta=5.65; %Параметр неравномерности функции Кайзера

fc1=125/FN; %Нормированные частоты среза

fc2=275/FN;

FC=[fc1 fc2]; % Вектор краевых частот

hn=fir1(N-1, FC, kaiser(N, beta)); %Получить коэффициенты фильтра

[H,f]=freqz(hn, 1, 512, FS);%Вычислить частотную характеристику

mag=20*log10(abs(H));

plot(f, mag), grid on

xlabel('Частота (Гц)')

ylabel('Амплитудная характеристика (дБ)')

Оптимизационные методы

Средство Signal Processing Toolbox в MATLAB содержит несколько программ разработки и функций для создания оптимальных КИХ–фильтров на основе алгоритмов Паркса–Мак–Клиллана и Ремеза. Основной командой для вычисления коэффициентов с помощью оптимального метода является remez. Команда может использоваться для разработки многополосных КИХ–фильтров с линейной фазовой характеристикой. В стандартной форме она имеет следующий синтаксис:

b = remez(N-l,F,M)

где N – длина фильтра, F – вектор нормированных граничных частот, а М – вектор желаемой амплитудной характеристики фильтра на заданных граничных частотах. Граничные частоты нормированы на половину частоты дискретизации и лежат в диапазоне от 0 до 1 (частота Найквиста соответствует 1).

Стандартную команду можно расширить, задав, например, относительные коэффициенты колебаний характеристики в полосе пропускания и подавления и/или тип фильтра. Если задаются относительные весовые коэффициенты и требуются амплитуды колебаний характеристики в полосе пропускания и полосе подавления, используется следующий синтаксис:

b = remez(N-l,F,M,WT),

где WT – вектор относительных коэффициентов амплитуд колебаний характеристики в полосах.

Для задания типа требуемого фильтра можно дополнительно добавить флаг 'ftype'. Существует четыре возможных типа, в зависимости от четности N и типа симметрии коэффициентов фильтра. Фильтры относятся к типу 1, когда длина фильтра нечетная (т.е. когда порядок фильтра, N–1, – четная величина), и к типу 2, когда длина фильтра четная. На использование фильтров типа 1 в качестве стандартных частотно-избирательных фильтров ограничений не существует. Фильтры типа 2 имеют нуль на частоте Найквиста, так что не могут использоваться при разработке фильтров верхних частот и заградительных фильтров. Фильтры типа 3 (N – нечетное) дают реализации преобразования Гильберта, и фильтры типа 4 (N – четное) - дифференциаторы. Для выбора фильтров типа 1 и 2 достаточно задать длину фильтра, но, чтобы указать тип фильтра для фильтров типа 3 и 4, необходимо включить флаг 'hilbert' или 'differentiator'.

Пример 3. Используя оптимальный метод, вычислите коэффициенты и изобразите частотную характеристику полосового КИХ–фильтра с линейной фазовой характеристикой, удовлетворяющего следующим условиям:

полоса пропускания 1000 – 1500 Гц

полоса перехода 500 Гц

длина фильтра 41

частота дискретизации 10 кГц

Решение

Полосы фильтра: 0–500 Гц (нижняя полоса подавления), 1000 –1500 Гц (полоса пропускания), 2000–5000 (верхняя полоса подавления). Граничные частоты нормируются на половину частоты дискретизации:

500/5000=0,1

1000/5000=0,2

1500/5000=0,3

2000/5000=0,4

5000/5000=1

Следовательно, вектор нормированных граничных частот F записывается в таком виде:

F=[0; 0,1; 0,2; 0,3; 0,4; 1]

Нужная амплитудная характеристика равна 1 в полосе пропускания и 0 в полосе подавления, поэтому вектор желаемой амплитудной характеристики записываем в виде

M=[001100]

Команды MATLAB:

Fs=10000; %Частота дискретизации

N=41; %Длина фильтра

M=[0 0 1 1 0 0]; %Желаемая амплитудная характеристика

F=[0 0.1 0.2 0.3 0.4 1]; %Края полос

b = remez(N-1, F, M); %Вычислить коэффициенты фильтра

[H, f] = freqz(b, 1, 512, Fs); %Вычислить частотную характеристику

mag = 20*log10(abs(H)); %фильтра и нарисовать ее

plot(f, mag), grid on

xlabel('Частота (Гц)')

ylabel('Амплитудная характеристика (дБ)')

Пример 4. Требуется полосовой фильтр с линейной фазовой характеристикой, удовлетворяющий следующим спецификациям:

полоса пропускания 12–16 кГц

ширина полосы перехода 2 кГц

неравномерность в полосе пропускания 1 дБ

затухание в полосе подавления 45 дБ

частота дискретизации 50 кГц

Оцените длину фильтра N и используйте оптимальный метод, чтобы определить его коэффициенты и амплитудно-частотную характеристику. Сравните амплитуды колебаний характеристики в полосе пропускания и подавления с заданными значениями.

Решение

Вначале определяем граничные частоты, затем нормируем их на половину частоты дискретизации:

10/25=0,4

12/25=0,48

16/25=0,64

18/25=0,72

Используя эти значения, получаем следующий вектор граничных частот:

F=[0, 0.4, 0.48, 0.64, 0.72, 1]

M=[0 0 1 1 0 0]

Для оценки длины фильтра можно использовать команду remezord. Для этого неравномерность в полосе пропускания и полосе подавления нужно перевести из децибелов в стандартные линейные единицы:

На основе граничных частот, желаемой амплитудной характеристики, значений амплитуды колебаний и частоты дискретизации оценивается порядок (N–1), а следовательно, и длина фильтра.

Команды MATLAB:

Fs=50000; %Частота дискретизации

Ap=1; %Неравномерность в полосе пропускания в дБ

As=45; %Затухание в полосе подавления

M=[0 1 0]; %Желаемая амплитудная характеристика

F=[10000, 12000, 16000, 18000] ; %Края полос

dp=(10^(Ap/20)-1)/(10^(Ap/20)+1); %Неравномерность в полосе пропускания и подавления

ds=10^(-As/20);

dev=[ds dp ds];

[N1, F0, M0, W] = remezord(F, M, dev, Fs);%Определить порядок фильтра

[b delta] = remez(N1, F0, M0, W); %Вычислить коэффициенты фильтра

[H, f] = freqz(b, 1, 1024, Fs); %Вычислить частотную характеристику

mag = 20*log10(abs(H)); %фильтра и нарисовать ее

plot(f, mag), grid on

xlabel('Frequency (Hz)')

ylabel('Magnitude (dB)')

Альтернативная программа. Порядок N=44. Используются нормированные частоты:

N=44; % Длина фильтра

Fs=50000; %Частота дискретизации

Ap=1; %Неравномерность в полосе пропускания в дБ

As=45; %Затухание в полосе подавления в дБ

M=[0 0 1 1 0 0]; %Желаемая амплитудная характеристика

F=[0, 0.4, 0.48, 0.64, 0.72 1]; % Края полос

dp=(10^(Ap/20)-1)/(10^(Ap/20)+1);

ds=10^(-As/20);

W=[dp/ds, 1, dp/ds];

dev=[ds ds dp dp ds ds ];

[b delta] = remez(N-1, F, M, W); %Вычислить коэффициенты фильтра

[H, f] = freqz(b, 1, 1024, Fs); %Вычислить частотную характеристику фильтра

mag = 20*log10(abs(H)); % и нарисовать ее

plot(f, mag), grid on

xlabel('Frequency (Hz)')

ylabel('Magnitude (dB)')

Метод частотной выборки

Для разработки КИХ–фильтров с произвольными частотными характеристиками, подобных фильтрам, рассмотренным при использовании метода частотной выборки, предназначена команда fir2. Синтаксис стандартной команды:

b=fir2(N-1, F,H)

Команда fir2 вычисляет коэффициенты КИХ–фильтра длины N. Вектор F задает нормированные частоты в интервале от 0 до 1 (причем частоты, как и ранее, нормированы на половину от частоты дискретизации). Вектор H определяет желаемую амплитудную характеристику в точках, заданных в F. Оба вектора должны иметь одинаковую длину.

Пример 5. КИХ–фильтр частотной выборки с линейной фазовой характеристикой имеет две частотные выборки в полосе перехода. Предполагается, что фильтр имеет 15 отводов и характеризуется следующими частотными выборками:

|H(k)|=1 k =0,1,2,3

|H(k)|=0,5571 k=4

|H(k)|=0,0841 k=5

|H(k)|=0 k=6,7

Определите коэффициенты фильтра, если частота дискретизации равна 2 кГц.

Решение

Частотные выборки заданы в диапазоне от 0 до половины частоты дискретизации. Следовательно, имеем такие положения точек выборки, нормированные на половину частоты дискретизации: 0, 1/7, 2/7, 3/7, 4/7, 5/7, 6/7/, 1.

Команды MATLAB:

Fs=2000; % Частота дискретизации

N=15; % Длина фильтра

fd=[0 1/7 2/7 3/7 4/7 5/7 6/7 1]; % Точки выборки

Hd=[1 1 1 1 0.5571 0.0841 0 0]; % Частотные выборки

hn=fir2(N-1, fd, Hd); % Вычислить импульсную характеристику

[H, f] = freqz(hn, 1, 512, Fs); % Расчет АЧХ и

plot(f, abs(H)), grid on

xlabel('Frequency (Hz)')

ylabel('Magnitude ')

ЗАДАНИЯ ДЛЯ ВЫПОЛНЕНИЯ

1. Выполните фильтрацию сигнала (табл. 1) двумя способами (с помощью conv и filter). Есть ли разница? Почему? Для расчета использовать КИХ-фильтр с коэффициентами: [0.0055 -0.0451 0.0692 -0.0554 -0.0634 0.5789 0.5789 -0.0634 -0.0554 0.0692 -0.0451 0.0055]

2. С помощью метода взвешивания требуется разработать 41 – точечный полосовой КИХ–фильтр, аппроксимирующий следующую идеальную амплитудную характеристику: H(f)=1, 2кГц ≤ f ≤ 4 кГц. На других частотах H(f) = 0. Частота дискретизации 10 кГц. Определите коэффициенты частотной характеристики фильтра и изобразите его частотные характеристики для следующих случаев:

2.1. прямоугольная весовая функция

2.2. весовая функция Хэмминга

3. Требуется, чтобы КИХ–фильтр нижних частот с линейной фазовой характеристикой удовлетворял следующим спецификациям: таблица 4. Используйте оптимальный метод.

4. Требуется, чтобы полосовой КИХ–фильтр удовлетворял следующим спецификациям: таблица 5. Используйте оптимальный метод.

5. Необходим полосовой КИХ–фильтр с линейной фазовой характеристикой, удовлетворяющий следующим спецификациям: таблица 6. Получите амплитудно-частотные характеристики фильтра для каждого из указанных случаев:

5.1. Используется весовая функция Хэмминга

5.2. Используется весовая функция Кайзера

5.3. Используется оптимальный метод

5.4. Сравните 3 варианта

6. КИХ-фильтр (N=15) частотной выборки с линейной фазовой характеристикой имеет три частотные выборки в полосе перехода и характеризуется следующими частотными выборками:

|H(k)|=1 k =0,1,2

|H(k)|=0,56 k=3

|H(k)|=0,21 k=4

|H(k)|=0,04 k=5

|H(k)|=0 k=6,7

Определите коэффициенты фильтра, если частота дискретизации равна 3 кГц.

7. Требуется многополосный КИХ–фильтр, удовлетворяющий следующим условиям: таблица 7. Используйте оптимальный метод для расчета коэффициентов фильтра и отобразите его АЧХ. Частоту дискретизации считайте равной 12 кГц, а ширину полосы перехода – 100 Гц.

КОНТРОЛЬНЫЕ ВОПРОСЫ

1. Как влияет на фазу сигнала фильтр с нелинейной фазовой характеристикой?

2. Напишите условия, при которых у фильтра будет линейная фазовая характеристика.

3. Определите типы фильтров с линейной фазовой характеристикой.

4. Фильтры какого типа не используются в качестве фильтров верхних частот? Почему?

5. Фильтры какого типа не используются в качестве фильтров нижних частот? Почему?

6. Фильтры какого типа наиболее универсальны?

7. Перечислите этапы разработки КИХ-фильтра.

8. Кратко опишите метод взвешивания.

9. Назовите несколько наиболее распространенных весовых функций.

10. Кратко опишите оптимизационный метод.

11. Кратко опишите метод частотной выборки.

12. Сравните метод взвешивания, оптимальный метод и метод частотной выборки.

Требования к содержанию отчета:

Составьте отчет в электронном виде, включающий все команды, вводимые в командной строке MATLAB и амплитудные характеристики фильтров.

Таблица 4

Фильтр нижних частот

Вариант	Длина фильтра	Край полосы пропускания, кГц	Край полосы подавления, кГц	Частота дискретизации, кГц
1	21	2	3	10
2	19	4	7	16
3	17	3	4	10
4	23	6	8	20
5	21	3	5	12
6	19	5	7	16
7	15	3	4	10
8	17	4	6	14
9	21	3	4	10
10	19	5	6	14
11	15	2	3	8
12	13	3	4	10
13	31	5	7	18
14	25	4	5	16
15	27	6	7	20

Таблица 5

Полосовой КИХ–фильтр

Вариант	Полоса пропускания, Гц	Ширина перехода, Гц	Неравномерность в полосе пропускания, дБ	Затухание в полосе подавления, дБ	Частота дискретизации, кГц
1	150-250	50	0.1	60	1
2	200-300	40	0.5	50	1
3	100-200	45	0.3	55	0.5
4	220-300	20	0.5	60	0.7
5	400-500	60	0.1	70	1.2
6	300-600	50	0.3	50	1.5
7	170-290	60	0.2	55	0.8
8	230-450	30	0.1	45	1
9	210-310	55	0.8	50	0.8
10	400-600	70	1	70	1.5
11	200-450	60	1	65	1.2
12	160-340	40	0.2	60	1
13	280-410	80	0.3	45	1.4
14	330-480	40	0.1	50	1.2
15	240-510	60	0.5	70	1.4

Таблица 6

Полосовой КИХ–фильтр

Вариант	Полоса пропускания, кГц	Неравномерность полосе подавления	Неравномерность в полосе пропускания	Частота дискретизации, кГц	Ширина перехода, кГц
1	8-12	0.001	0.01	48	3
2	7-13	0.002	0.05	40	2
3	10-20	0.003	0.03	64	3
4	13-20	0.004	0.05	48	2
5	10-14	0.005	0.01	32	1
6	8-10	0.001	0.03	24	1
7	6-12	0.002	0.02	44	3
8	3-6	0.003	0.01	20	1
9	12-18	0.004	0.08	64	2
10	8-20	0.005	0.01	50	3
11	7-11	0.001	0.01	26	3
12	12-22	0.002	0.02	64	2
13	10-13	0.003	0.03	44	1
14	11-18	0.004	0.01	44	2
15	3-5	0.005	0.05	16	3

Таблица 7

Многополосный КИХ–фильтр

Вар	Полоса 1, кГц	A_s1 дБ	A_p₁ дБ	Полоса 2, кГц	A_s2 дБ	A_p2 дБ	Полоса 3, кГц	A_s3 дБ	A_p3 дБ	Полоса 4, кГц	A_s4 дБ	A_p4 дБ	Полоса 5, кГц	A_s₅ дБ	A_p5 дБ
1	0 – 0.5	49	0.2	1-1.5	38	0.3	1.8– 2.5	38	0.2	3 – 3.6	35	0.3	4.1–5	55	0.1
2	0-0.3	45	0.3	1.1-1.3	45	0.1	1.7-2.1	45	0.3	2.9-3.4	40	0.2	4.5-5	60	0.2
3	0-0.4	47	0.2	1.2-1.4	40	0.2	1.8-2.4	40	0.2	2.8-3.5	45	0.5	4.2-5	40	0.4
4	0-0.6	38	0.1	0.9-1.5	50	0.5	1.9-2.6	50	0.1	3-3.7	48	0.1	4.1-5	50	0.1
5	0-0.5	40	0.3	0.9-1.5	35	0.3	1.8-2.4	35	0.3	2.9-3.4	30	0.4	4.2-5	60	0.3
6	0-0.3	50	0.3	1-1,5	47	0.2	1.8-2.1	47	0.3	3.1-3.8	44	0.7	4-5	44	0.5
7	0-0.4	30	0.5	1-1.4	35	0.7	1.9-2.5	35	0.5	2.9-3.4	33	0.5	4.1–5	32	0.5
8	0-0.3	40	0.2	1.1-1.3	36	0.8	1.7-2.1	36	0.2	3.1-3.8	40	0.4	4.5-5	58	0.5
9	0-0.6	50	0.1	1-1,5	46	0.1	1.8-2.1	46	0.1	2.8-3.5	43	0.1	4.2-5	32	0.2
10	0-0.5	44	0.7	0.9-1.5	36	0.5	1.9-2.5	36	0.7	3-3.7	34	0.9	4-5	57	0.7
11	0-0.3	46	0.1	1.2-1.4	41	0.4	1.8-2.1	41	0.1	2.9-3.4	43	0.2	4.3 –5	54	0.2
12	0-0.4	39	0.3	1-1,5	53	0.3	1.9-2.5	53	0.3	3.1-3.8	50	0.8	4.1-5	45	0.7
13	0-0.6	51	0.2	1.2-1.4	60	0.3	1.7-2.1	60	0.2	2.8-3.5	60	0.3	4.2-5	51	0.3
14	0-0.5	57	0.1	1.1-1.3	55	0.7	1.9-2.7	55	0.1	3-3.7	55	0.1	4.5-5	33	0.2
15	0-0.3	39	0.1	1 - 1,5	43	0.1	1.8-2.1	43	0.1	2.9-3.4	41	0.2	4.1–5	37	0.5

A_p - неравномерность в полосе пропускания

A_s - Затухание в полосе подавления

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Айфичер Э.С., Джервис Б.У. Цифровая обработка сигналов: практический подход [Текст]: Пер. с анг. – М.: Издательский дом «Вильямс», 2004. – 992с.

2. Сергиенко А.Б. Цифровая обработка сигналов [Текст] – СПб.: Питер, 2002. – 608с.

Дата добавления: 2018-04-05; просмотров: 4834; Мы поможем в написании вашей работы!

Поделиться с друзьями:

⇐ Предыдущая 1 2 3 4 5 6 7 89

Мы поможем в написании ваших работ!