CC BY
229Моделирование цифровом обработки сигналов в MATLAB
Часть 9. Моделирование цифровых преобразователей Гильберта и дифференциаторов программными средствами и средствами GUI FDATool MATLAB
В цикле «Моделирование цифровой обработки сигналов в MATLAB» предыдущие статьи [6—13] были посвящены моделированию цифровых фильтров (ЦФ), КИХ и БИХ, в том числе с фиксированной точкой (ФТ), программными средствами и средствами GUI (Graphic User Interface — графический интерфейс пользователя) MATLAB .
Настоящая статья завершает тему моделирования цифровых фильтров и посвящена цифровым преобразователям Гильберта и дифференциаторам, которые в цифровой обработке сигналов реализуются на основе КИХ-фильтров с линейной ФЧХ (ЛФЧХ). Рассматривается их синтез программными средствами и средствами GUI FDATool (Filter Design and Analysis Toolbox — средства проектирования и анализа фильтров) MATLAB.
Алла СОЛОНИНА
Цифровой преобразователь Гильберта
Цифровым преобразователем Гильберта (ЦПГ) называют линейную дискретную систему, формирующую на выходе пару дискретных сигналов, сопряженных по Гильберту (фазы сигналов отличаются на п/2) в заданной рабочей полосе [5].
ЦПГ может быть реализован на базе КИХ-фильтров 3-го и 4-го типов (табл. 1 в [6]), ЛФЧХ которых обеспечивает сдвиг фазы на п/2. Предпочтение отдается КИХ-фильтру 3-го типа, так как он позволяет получить импульсную характеристику (ИХ) к(н), каждый второй отсчет которой равен нулю, тем самым сокращается число арифметических операций при вычислении реакции ЦПГ, что весьма важно при его реализации, например, на цифровом процессоре обработки сигналов (ЦПОС).
На базе КИХ-фильтра 3-го типа можно синтезировать только полосовой фильтр (ПФ),
при этом специфика требований к АЧХ ЦПГ, по сравнению с требованиями к АЧХ ПФ, будет следующей:
• АЧХ ЦПГ должна быть симметричной относительно середины ^/4 основной полосы частот [0; /д/2] [6] для получения ИХ к(н), каждый второй отсчет которой равен нулю. Поэтому требования к АЧХ ЦПГ задаются симметричными относительно _/д/4.
• Рабочая полоса ЦПГ не должна превосходить полосу пропускания ПФ.
• Максимально допустимое отклонение в рабочей полосе 8раб не должно быть меньше максимально допустимого отклонения 81 в ПП (как правило, с ним совпадает 8раб = 81).
• Максимально допустимое отклонение 82 в ПЗ нет необходимости задавать слишком жестко, так как эффективность ЦПГ оценивается в рабочей области.
Пример 1
Синтезировать ЦПГ на базе ПФ КИХ-фильтра 3-го типа методом чебышевской ап-
проксимации [6] с автоматическим подбором минимального порядка при следующих требованиях:
• частота дискретизации f = 1000 Гц;
• левая граничная частота рабочей полосы
граблен = 100 Гц;
• правая граничная частота рабочей полосы
/раб_прав = 400 Гц;
• максимально допустимое отклонение в рабочей полосе 8раб = Si = 0,05;
• максимально допустимое отклонение в ПЗ1 и ПЗ2 S2 = 0,1.
Требования к АЧХ ПФ, на базе которого синтезируется ЦПГ, приведены в таблице.
По требованиям к АЧХ будем синтезировать ЦПГ (ПФ) минимального порядка с помощью функции firgr [6] на базе КИХ-филь-тра 3-го типа ('hilbert') с параметром m, равным 'mineven' (аналогичный пример 5 в [6]):
>> Fs=1000;
>> fki=50; fti=100; ft2=400; fk2=450; f=[fk1 ft1 ft2 fk2];
>> m=[0 1 0];
>> d2=0.1; d1=0.05; ripple=[d2 d1 d2];
>> [R,f0,m0,weight]=firpmord(f,m,ripple,Fs);
>> f0' ans =
0 0.1000 0.2000 0.8000 0.9000 1.0000
>> m0' ans =
0 0 1 1 0 0
>> weight' ans =
1 2 1
>> [b,err,opt]=firgr({'mineven',R},f0,m0,ripple,'hilbert');
>> opt.order
ans =
22
Выведенные значения f0, m0 и weight будут использованы далее в примере 4.
Таблица. Требования к АЧХ ПФ при синтезе ЦПГ
Частоты (Гц) и их обозначения в MATLAB Максимально допустимые отклонения АЧX и их обозначение в MATLAB
Частота дискретизации Fs 1000 - - -
Граничная частота ПЗ1 f-k fk1 50 Максимально допустимое отклонение от нуля в ПЗ1 82 d2 0,1
Граничная частота ПП f-L ft1 100 Максимально допустимое отклонение от единицы в ПП Si d1 0,05
Граничная частота ПП fX ft2 400
Граничная частота ПЗ2 fk fk2 450 Максимально допустимое отклонение от нуля в ПЗ2 82 d2 0,1
Синтезирован оптимальный ПФ с ЛФЧХ порядка Л ( = 22 на базе КИХ-фильтра 3-го типа, который функционирует как ЦПГ в рабочей полосе.
Построим графики ИХ, АЧХ и ФЧХ ЦПГ с помощью внешней функции р1о^к [6] (рис. 1):
Каждый второй отсчет ИХ ЦПГ равен нулю. На рис. 2 представлена структура ЦПГ [5]. Рассмотрим работу ЦПГ при входном гармоническом сигнале:
x (n) = cos((2nf/^)n),
где f = 125 Гц. Выходной сигнал y(n) рассчитаем с помощью функции filter [11], для которой вектор коэффициентов b был рассчитан в примере 1 (рис. 3):
>> Fs=1000; f=125; n=0:63;
>> x=cos((2.*pi.*f./Fs).*n);
>> subplot(2,1,1), stem(n,x,'fill','MarkerSize',3),... xlabel('n'), title('x(n)'), grid >> a=[1];
>> y=filter(b,a,x);
>> subplot(2,1,2), stem(n,y,'fin','MarkerSize',3),... xlabel('n'), title('y(n)'), grid
Согласно рис. 2, при N = 23 сопряженными по Гильберту должны быть сигналы x[n-(N-1)/2] = x(n-11) и y(n). Для формирования задержанного сигнала x(n-11) (далее обозначен как xx) создадим внешнюю функцию shift [5]:
function [y,k]=shift(x,n,m)
% Сдвиг последовательности % n — дискретное время n=0,1,2...
% x — не задержанная последовательность x(n)
% m — задержка % k — дискретное время k=n+m
% у — задержанная последовательность y(k)=x(n-m);
k=n+m;
y=x;
Графические результаты представлены на рис. 4:
>> [xx,k]=shift(x,n,11);
>> subplot(2,1,1), stem(k,xx,'fill','MarkerSize',3),... xlabel('n'), title('x(n-11)'), grid >> subplot(2,1,2), stem(n,y,'fin','MarkerSize',3),... xlabel('n'), title('y(n)')
Наконец, для того чтобы убедиться в сопряжении по Гильберту сигналов x (n-11) и y(n), сравним их в установившемся режиме на одинаковом интервале, например, n = [27; 63] (рис. 5):
>> subplot(2,1,1), stem(k,xx,'fill','MarkerSize',3),... xlabel('n'), title('x(n-11) n=[27;63]'), xlim([27,63]), grid >> subplot(2,1,2), stem(n,y,'fin','MarkerSize',3),... xlabel('n'), ylabel('y(n) n=[27;63]'), xlim([27,63]), grid
Как видим, фазы входного и выходного гармонических сигналов отличаются на п/2.
>> R=opt.order; >> plot_fir(R,b,Fs)
н
0,5
-0,5
x(n—11) n = [27; 63]
ЗО 35 40 45 50 55 60
Рис. 5. Сигналы ЦПГ в установившемся режиме: а) входной задержанный; б) выходной
Цифровой дифференциатор
Цифровым дифференциатором (ЦД) называют линейную дискретную систему, обеспечивающую вычисление отсчетов производной сигнала по отсчетам самого сигнала в рабочей полосе, совпадающей с основной полосой частот либо расположенной внутри нее.
Идеальный ЦД имеет в рабочей полосе частотную характеристику:
_ “і?
Я(еу“) = ]козе] 2 =
= к&е 12 2>=А{іи)еЖш)
с ЛФЧХ, соответствующей КИХ-фильтрам 3-го и 4-го типов (табл. 1 в [6]).
Идеальная АЧХ ЦД в рабочей полосе в шкале абсолютных частот имеет вид:
А(/) = k (//(/д/2)Х при
Реальный ЦД синтезируется на основе КИХ-фильтров 3-го или 4-го типа по требованиям к АЧХ в рабочей полосе [_/раб1; /рабг]> совпадающей с основной полосой [0; /д/2] либо расположенной внутри нее.
В зависимости от расположения рабочей полосы в основной полосе [0; /д/2] различают следующие разновидности ЦД [5]:
• широкополосные ЦД (ШЦД) — с рабочей полосой [0; /д/2];
• низкочастотные ЦД (НЦД) — с рабочей полосой [0; /раб];
• полосовые ЦД (ПЦД) — с рабочей полосой [/раб1; /раб2];
• высокочастотные ЦД (ВЦД) — с рабочей
полосой [/раб; /д/2].
Так как АЧХ КИХ-фильтров 3-го типа по определению равна нулю на границах основной полосы частот, ШЦД и ВЦД можно синтезировать только на базе КИХ-фильтра 4-го типа.
Максимально допустимое отклонение АЧХ от идеальной в рабочей полосе (8раб) задается достаточно жестко, а в нерабочей (8нераб) — более или менее формально, 8нераб^8раб.
ЦД синтезируются с помощью функций firpm или firgr [6] с одинаковыми форматами:
[b,error,opt]=firpm(R,f0,m0,weight,'differentiator',{lgrid})
[b,error,opt]=firgr(R,f0,m0,weight,'differentiator',{lgrid})
Здесь R — оценка порядка фильтра R (формула (2) в [6]), задаваемая произвольно, но с учетом типа КИХ-фильтра — 3-го или 4-го; f0 — вектор нормированных частот f = //(/д/2) в основной полосе частот [0; 1]; длина вектора f0 — всегда четное число, включающее следующие частоты в порядке их следования слева направо:
• левую границу основной полосы частот f0 = 0;
• граничные частоты рабочей полосы: отсутствуют — для ШЦД; одна—для НЦД и ВЦД; две — для ПЦД;
• граничные частоты нерабочей полосы: отсутствуют — для ШЦД; одна — для НЦД и ВЦД; две — для ПЦД;
• правую границу основной полосы частот /0 = 1;
m0 — вектор значений идеальной АЧХ ЦД на частотах вектора f0; длины векторов f0 и m0 совпадают. При задании значений вектора m0 необходимо помнить о следующем: m0 = 0 на частоте f0 = 0; m0 = 0 на частоте f0 = 1 для КИХ-фильтров 3-го типа. weight — вектор весов в рабочей и нерабочей полосах в порядке их следования слева направо. Длина вектора weight вдвое меньше длины вектора f0 и равна: 1 — для ШЦД; 2 — для НЦД и ВЦД; 3 — для ПЦД.
Веса рассчитываются по методике, рассмотренной в [6], на основе заданных максимально допустимых отклонений 8раб и 8нераб. В рабочей полосе вес всегда равен 1, а в нерабочей 8нераб/8раб>1.
Примечания. В MATLAB граница нерабочей полосы обязательно задается и для ВЦД, и для ПЦД. Длина вектора f0 будет равна: 2 — для ШЦД; 4 — для НЦД и ВЦД; 6 — для ПЦД.
Параметры lgrid, error и opt для функции firpm были определены ранее в [6], но с одним существенным замечанием. Для ЦД модуль максимального отклонения реальной АЧХ от идеальной в рабочей полосе равен не значению параметра error, как для частотноизбирательных КИХ-фильтров, а значению параметра max(abs(opt.error)) (аналогично примеру 4 [6]).
Пример 2
Синтезировать ШЦД с частотой дискретизации /д = 8000 Гц при максимально допустимом отклонении АЧХ от идеальной в рабочей полосе 8раб = 0,02 и k = 1.
ШЦД будем синтезировать на базе КИХ-фильтра 4-го типа:
>> f0=[0 1];
>> m0=[0 1];
>> weight=[1];
>> R=15;
>> [b,error,opt]=firpm(R,f0,m0,weight,'differentiator');
Проверим выполнение требований к АЧХ ШЦД в рабочей полосе:
>> max(abs(opt.error)) ans =
0.0136
Требования выполняются. Последовательно уменьшая нечетный порядок фильтра, получим ШЦД минимального порядка Rmin =11.
Построим графики импульсной характеристики, АЧХ и ФЧХ ШЦД с помощью внешней функции plot_fir [6] (рис. 6):
>> Fs=8000;
>> plot_fir(R,b,Fs)
Рис. 8. АЧХ цифрового преобразователя Гильберта
Пример 3
Синтезировать НЦД с частотой дискретизации ^ = 8000 Гц при границе рабочей полосы _/раб = 2000 Гц £—), границе нерабочей частоты /нераб = 3000 Гц (£п—), максимально допустимом отклонении АЧХ от идеальной в рабочей полосе 8раб = 0,02, в нерабочей —
8нераб = 0,1 и к = 0,5
НЦД можно синтезировать на базе КИХ-фильтров 3-го и 4-го типов. Будем его синтезировать на базе КИХ-фильтра 3-го типа:
>> Fs=8000;
>> £-=2000; £п—=3000; f0=[0 fw/(Fs/2) fnw/(Fs/2) 1]
£0 =
0 0.5000 0.7500 1.0000
>> т0=[0 0.25 0.375 0];
>> ше1дЬ1=[5 1];
>> Я=12;
>> [b,error,opt]=£irpm(R,£0,m0,weight,'di££erentiator');
Проверим выполнение требований к АЧХ НЦД в рабочей полосе:
>> max(abs(opt.error))
ап =
0.0276
Требования не выполняются. Последовательно увеличивая порядок фильтра, получим НЦД минимального порядка Ктп = 14.
Построим графики ИХ, АЧХ и ФЧХ НЦД с помощью внешней функции рЬ^Пг [6] (рис. 7):
Синтез цифровых преобразователей Гильберта средствами GUI FDATool
Синтез ЦПГ средствами GUI FDATool предполагает обязательное знакомство с синтезом ЦПГ программными средствами MATLAB, рассмотренными выше.
Пример 4
В примере 1 был синтезирован ЦПГ программными средствами на базе КИХ-фильтра ПФ 3-го типа методом чебышевской аппроксимации при заданных требованиях к ЦПГ, на основе которых были сформулированы требования к АЧХ КИХ-фильтра ПФ (таблица). Синтезировать тот же ЦПГ средствами GUI FDATool.
Основные этапы синтеза ЦПГ включают в себя (рис. 8):
1. Выбор типа ЦФ. В группе Design Method — переключатель FIR.
2. Выбор метода синтеза. В группе Design Method в раскрывающемся списке FIR — Equiripple.
3. В раскрывающемся списке специальных фильтров (без имени) выбор Hilbert Transformer («Преобразователь Гильберта»).
4. Задание входных параметров [12] ЦПГ. Входные параметры ЦПГ рассчитываются программными средствами MATLAB с помощью функции firpmord по требованиям к АЧХ ПФ (КИХ-фильтру 3-го типа), на основе которого синтезируется ЦПГ. В данном случае требования к АЧХ ЦФ заданы в таблице, а входные параметры ЦПГ рассчитаны с помощью функции firpmord в примере 1. В GUI FDATool эти входные параметры задаются следующим образом:
• В группе Frequency and Magnitude Specifications:
- в раскрывающемся списке Frequency Units — Normalized (0 to 1) (по умолчанию);
- в поле ввода Freq. vector («Вектор частот») — [0 0.1 0.2 0.8 0.9 1];
- в поле ввода Mag. vector («Вектор АЧХ») — [0 0 1 1 00];
- в поле ввода Weight vector («Вектор весов») — [1 2 1].
Примечание. Векторы Freq. vector, Mag. vector и Weight vector тождественны соответственно векторам f0, m0 и weightв функции firpmord ирассчитаны с ее помощью.
• В группе Filter Order в списке ввода Specify order — порядок фильтра 18. Примечание. Порядок фильтра тождественен параметру R в функции firpmord и рассчитан с ее помощью.
• В группе Options в поле ввода Density Factor — 20.
Примечание. Параметр Density Factor тождественен параметру Igrid в функциях firpm и firgr [6].
5. Синтез ЦПГ производится после нажатия кнопки Design Filter. По завершении синтеза автоматически выдаются:
• В группе Magnitude Response (dB) — график АЧХ (дБ) — характеристика ослабления (5) в [6].
• В группе Current Filter Information:
- Structure — Direct-Form FIR;
- Order — 18;
- Stable — Yes;
- Source — Designed.
Используя соответствующие команды в пункте меню Analysis, легко убедиться, что у синтезированного ЦПГ:
• АЧХ симметрична относительно частоты /д/4 (нормализованной частоты 0,5).
• импульсная характеристика антисимметрична — это КИХ-фильтр 3-го типа — и каждый ее второй отсчет равен нулю. Отметим, что синтезированный ЦПГ порядка R = 18 не является оптимальным, в чем легко убедиться, проверив выполнение требований к АЧХ в ПП и ПЗ. Оптимальный ЦПГ будет синтезирован при Rmin = 22 (рис. 8).
Синтез цифровых дифференциаторов средствами GUI FDATool
Синтез ЦД средствами GUI FDATool предполагает обязательное знакомство с синтезом ЦД программными средствами MATLAB, рассмотренными выше.
>> plot_fir(R,b,Fs)
Рис. 9. АЧХ низкочастотного цифрового дифференциатора
Пример 5
В примере 3 был синтезирован низкочастотный цифровой дифференциатор (НЦД) программными средствами на базе КИХ-фильтра 3-го типа методом чебышевской аппроксимации. Синтезировать тот же НЦД средствами GUI FDATool.
Основные этапы синтеза НЦД включают в себя (рис. 9):
1. Выбор типа ЦФ. В группе Design Method — переключатель FIR.
2. Выбор метода синтеза. В группе Design Method в раскрывающемся списке FIR — Equiripple.
3. В раскрывающемся списке специальных фильтров (без имени) выбор Differentiator («Дифференциатор»).
4. Задание входных параметров [12] ЦД. Входные параметры НЦД тождественны входным параметрам функций firpm (или firgr), и они задаются в GUI FDATool следующим образом:
• В группе Frequency and Magnitude Specifications:
- в раскрывающемся списке Frequency Units — Normalized (0 to 1) (по умолчанию);
- в поле ввода Freq. vector («Вектор частот») — [0 0.5 0.75 1];
- в поле ввода Mag. vector («Вектор АЧХ») — [0 0.25 0.5 0];
- в поле ввода Weight vector («Вектор весов») — [5 1].
Примечание. Векторы Freq. vector, Mag. vector и Weight vector тождественны соответственно векторам f0, m0 и weight, рассчитанным с помощью функции firpm (пример 3);
• В группе Filter Order в списке ввода Specify order — порядок фильтра 12. Примечание. Порядок фильтра тождественен параметру R (пример 3);
• В группе Options в поле ввода Density Factor — 20.
Примечание. Параметр Density Factor тождественен параметру Igrid в функциях firpm и firgr [6].
5. Синтез НЦД производится после нажатия кнопки Design Filter. По завершении синтеза автоматически выдаются:
• В группе Magnitude Response (dB) — график АЧХ (дБ) — характеристика ослабления (5) в [6].
• В группе Current Filter Information:
- Structure — Direct-Form FIR;
- Order — 12;
- Stable — Yes;
- Source — Designed.
Отметим, что синтезированный НЦД порядка R = 14 не является оптимальным, в чем легко убедиться, проверив выполнение требований к АЧХ. Оптимальный НЦД будет синтезирован при Rmin = 14 (рис. 9).
Литература
1. Ingle V., Proakis J. Digital Signal Processing Using MATLAB. Second Edition. Thomson.
2. Оппенгейм А., Шафер Р. Цифровая обработка сигналов. М.: Техносфера, 2006.
3. Сергиенко А. Б. Цифровая обработка сигналов, 2-е изд. СПб.: ПИТЕР, 2006.
4. Солонина А. И., Улахович Д. А., Арбузов С. М., Соловьева Е. Б. Основы цифровой обработки сигналов. 2-е изд. СПб.: БХВ-Петербург, 2005.
5. Солонина А. И., Арбузов С. М. Цифровая обработка сигналов. Моделирование в MATLAB. СПб.: БХВ-Петербург, 2008.
6. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 1. Синтез оптимальных (по Чебышеву) КИХ-фильтров программными средствами MATLAB // Компоненты и технологии. 2008. № 11.
7. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 2. Синтез оптимальных БИХ-фильтров программными средствами MATLAB // Компоненты и технологии. 2008. № 12.
8. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 3. Описание структур КИХ- и БИХ-фильтров в MATLAB // Компоненты и технологии. 2009. № 1.
9. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 4. Моделирование структур цифровых фильтров с фиксированной точкой программными средствами MATLAB: анализ характеристик КИХ-фильт-ров // Компоненты и технологии. 2009. № 2.
10. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 5. Моделирование структур цифровых фильтров с фиксированной точкой программными средствами MATLAB: анализ характеристик БИХ-фильт-ров // Компоненты и технологии. 2009. № 3.
11. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 6. Моделирование структур цифровых фильтров с фиксированной точкой программными средствами MATLAB: квантование воздействия и вычисление реакции // Компоненты и технологии. 2009. № 4.
12. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 7. Моделирование цифровых фильтров средствами программ GUI MATLAB: GUI FDATool // Компоненты и технологии. 2009. № 5.
13. Солонина А. Моделирование цифровой обработки сигналов в MATLAB. Часть 8. Моделирование цифровой фильтрации средствами программ GUI MATLAB: GUI SPTool // Компоненты и технологии. 2009. № 6.
