Описать методику проектирования КИХ-фильтра с помощью окон.
Выбрать критерии точности фильтрации детерминированного и случайного сигнала КИХ-фильтром.
Составить схему определения точности фильтрации в приложении Simulink / Профильтровать разные виды сигналов тремя фильтрами: выбранным для взвешивания, спроектированным с помощью прямоугольного и спроектированным с помощью хорошего окна в MATLAB.
По результатам моделирования найти зависимости точности фильтрации для разных входных сигналов, уровня боковых лепестков, типов окон и т. п..

Методика проектирования КИХ-фильтра с помощью окон

1. Понятие об эквивалентности дискретной и непрерывной систем

Пусть имеется непрерывная система с передаточной функцией и импульсной характеристикой . Они связаны преобразованием Фурье. Пусть имеется дискретная система с передаточной функцией и импульсной характеристикой .

Дискретная система считается эквивалентной системе непрерывной, если в отсчётах соблюдается равенство сигналов на выходе систем при одинаковых входных данных.

Входные данные:

(Формула 1.1)

Для эквивалентности двух систем необходимо и достаточно, чтобы импульсная характеристика непрерывной системы равнялась импульсной характеристике дискретной системы в конкретные отсчёты времени:

(формула 1.2)

В пакете MATLAB для определения коэффициентов импульсной передаточной функции из непрерывного аналога служит функция

[bz,az] = impinvar(b,a,fs)

Где bz, az – вектора коэффициентов числителя и знаменателя импульсной передаточной функции , b, a – вектора коэффициентов числителя и знаменателя передаточной функции непрерывной системы , а fs – значение частоты дискретизации [Гц].

В качестве примера действия функции рассмотрим конвертирование аналогового фильтра нижних частот в цифровой фильтр с частотой дискретизации в 10 Гц:

[b,a] = butter(4,0.3,'s');
[bz,az] = impinvar(b,a,10)
bz = 1.0e-006 * -0.0000 0.1324 0.5192 0.1273 0
az =1.0000 -3.9216 5.7679 -3.7709 0.9246

2. Умножение импульсной характеристики на прямоугольное окно

Проектирование цифровых фильтров с помощью окон основывается на понятии об эквивалентной импульсной характеристике.

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

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

Однако для физических систем импульсная характеристика неограниченна во времени, что не позволяет так просто решить задачу проектирования фильтра. Чтобы не проектировать физически нереализуемый фильтр, с обоих сторон «обрезают» число отсчётов импульсной характеристики фильтра. Это эквивалентно умножению импульсной характеристики фильтра на окно, простейшее из которых – прямоугольное ().

(формула 1.3)

- фактическая импульсная характеристика, с помощью которой можно проектировать физически реализуемый фильтр.

Импульсная характеристика цифрового фильтра: перемножение импульсной характеристики идеального фильтра на импульсную характеристику прямоугольного окна.

Рис. 1.1. Умножение импульсных характеристик дискретного фильтра и окна.

Частотная характеристика синтезируемого фильтра есть свертка частотной характеристики идеального фильтра с частотной характеристикой окна.

Спектр произведения равен свёртке спектров множителей:

(формула 1.4)

Для оценки влияния прямоугольного окна находят частотную передаточную функцию окна:

(формула 1.5)

В пакете MATLAB прямоугольное окно моделируется функцией w = rectwin(L)

Где L – ширина прямоугольного окна являющаяся также порядком фильтра,

w – вектор-столбец значений окна (вычисляются значения ординаты через две временных еденицы L).

У этого метода есть свои существенные недостатки:

Передаточная функция прямоугольного окна имеет боковые лепестки, уровень которых достигает 20%. Точный характер искажения определяется при задании конкретной весовой функции и вычислении свёртки (формула 1.5). Очевидно, что спектр результата будет размываться.
Спектр окна бесконечен, но при попытке сильно его ограничить возникает эффект Гиббса – волнистость и резкие всплески. Обычно подобные проблемы легко устраняются увеличением числа гармоник, описывающих функцию, но на данный эффект этот метод почти не действует, хотя колебания становятся более высокочастотными, а выбросы более короткими. Для уменьшения эффекта Гиббса используют окна с плавно спадающими до нуля краями. Для их описания требуется меньшее количество гармоник, а влияние боковых лепестков становится менее заметным.

3. Эффект Гиббса

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

Знакомство с явлением Гиббса мы начнем с рассмотрения простейшей разрывной функции

(формула 1.6)

Известно, что для всех справедливо равенство

(формула 1.7)

Можно показать, что график функции

(формула 1.8)

на промежутке колеблется вокруг значения у=1, причём амплитуда горбов убывает по мере удаления от концов промежутка, а глобального максимума функция S₂_n_-1достигает в точках и .

С учётом сказанного, для амплитуды функции S₂_n_-1на промежутке будем иметь

(формула 1.9)

В последнем выражении легко узнать сумму Римана интеграла:

(формула 1.10)

соответствующую разбиению интервала на n равных частей.

(формула 1.11)

Входит в значение амплитуды при . По таблицам находим

Амплитуда частотной суммы S_2n-1 ряда Фурье функции f примерно на 18% превышает амплитуду самой функции f.

Более тонкий анализ показывает, что эффект Гиббса имеет место (причем с той же самой постоянной 18%) для любой кусочно-гладкой функции f.

4. Window Design & Analysis Tool

Окна характеризуются своей формой и амплитудно-частотной характеристикой. Высокий уровень амплитуды (называемый главным лепестком) характеризует полосу пропускания, низкий уровень амплитуды (называемый боковыми лепестками) характеризует полосу подавления. Между главным лепестком и боковыми лепестками существует переходная область. Чем она уже, тем лучше подавляются ненужные частоты сигнала, но, вместе с тем, увеличивается порядок фильтра, который также называют длительностью импульсной характеристики (параметр L). Порядок фильтра определяет количество арифметических операций над отсчетами и объем сигнальной памяти при реализации фильтра. Чем длиннее будет импульсная характеристика, тем большая нагрузка будет ложиться на аппаратные ресурсы.

Проектирование фильтра с конечной импульсной характеристикой в пакете прикладных программ для решения задач технических вычислений MATLAB с помощью окон осуществляется с помощью программы синтеза и анализа весовых (оконных) функций Window Design and Analysis Tool.

Программа вызывается посредством команды

>> wintool

И является графической оболочкой (GUI, Graphical user interface, Графический Интерфейс Пользователя) для вызова имеющихся в пакете функций расчёта окон.

Signal Processing Toolbox: Window Design and Analysis Tool

Рисунок 1.2. Окно GUI Window Design and Analysis Tool.

На окне Window Design and Analysis Tool (Рисунок 1.2) имеется три панели:

· На первой (верхней) панели Window Viever показаны временные отсчёты и отсчёты по частоте. Отсчёты отображаются в соответствии с данными, введёнными в двух нижних панелях. В один момент времени могут быть показаны отсчёты только одного окна (на рисунке 1.2 окно носит имя window_1), выбранного из всех созданных пользователем в левой нижней панели Window List, список окон. Ниже двух графиков показаны три величины, связанные с окнами:

· Leakage Factor, фактор утечки – отношение мощности боковых лепестков к полной мощности окна, указывается в процентах.

· Relative sidelobe attenuation, относительное ослабление лепестков – различие по высоте первого пика [dB] (в английской литературе – mainlobe, «главный лепесток») и наибольшего бокового лепестка [dB].

· Mainlobe width (-3dB), ширина главного лепестка, уменьшенная на 3 дБ. (Дело в том, что границу полосы пропускания фильтра [dB] в проектировании с помощью окон принимают равной наибольшей ординате второго графика минус три дБ, пояснение на рисунке 1.3)

Signal Processing Toolbox: Window Design and Analysis Tool - Frequency domain

Рис. 1.3. Граница полосы пропускания в окне Window Design and Analysis Tool.

· На второй панели Window List можно создавать и выбирать (из уже созданных) окна (Рисунок 1.4). Команда Add a new window добавляет новое окно Хемминга (по умолчанию) с параметром длины по временной оси в 64 единицы и именем window_КоличествоУжеСозданныхОкон+1 в список. Команда Copy window копирует окно window_n в copy_of_window_n, этот пункт удобен для сравнения окон. Команда Delete удаляет выделенное окно из списка. Команда Save to workspace сохраняет отмеченные окна в виде вектора (сохраняются значения Amplitude графика Time domain через две временных еденицы), причём имя окна в wintool совпадает с именем полученного вектора (далее с помощью команды window_1 продемонстрирован вектор значений для окна window_1 рисунка 1.2). Если в списке, удерживая Ctrl или Shift, выделить несколько проектируемых окон, то на графиках они отобразятся разными цветами, что очень удобно для их сравнения.

Signal Processing Toolbox: Window Design and Analysis Tool - Window List

Рис. 1.4. Список окон Design and Analysis Tool.

>> window_1 =

0.0004 0.0009 0.0018 0.0034 0.0060 0.0097 0.0150 0.0224 0.0323 0.0452 0.0617 0.0821 0.1071 0.1368 0.1717 0.2119 0.2573 0.3076 0.3626 0.4214 0.4834 0.5473 0.6121 0.6763 0.7386 0.7972 0.8509 0.8982 0.9377 0.9684 0.9894 1.0000 1.0000 0.9894 0.9684 0.9377 0.8982 0.8509 0.7972 0.7386 0.6763 0.6121 0.5473 0.4834 0.4214 0.3626 0.3076 0.2573 0.2119 0.1717 0.1368 0.1071 0.0821 0.0617 0.0452 0.0323 0.0224 0.0150 0.0097 0.0060 0.0034 0.0018 0.0009 0.0004

В целях экономии места из текста убраны переводы на новую строку.

· На третьей панели, Current Window Information (Рисунок 1.5), текущая информация окна, указываются параметры окна, такие как Name (имя), Type (тип фильтра), MATLAB Code (код для окна, созданного пользователем, User Defined), Length (временной интервал), Parameter (специальный параметр для фильтра: SidelobeAtten для фильтра Чебышева, Alpha для фильтров Гаусса и Тьюки, Beta для фильтра Кайзера), Sampling (указывается для фильтров с косинусоидальным сглаживанием).

Signal Processing Toolbox: Window Design and Analysis Tool - Current Window Information

Рис. 1.5. Панель Current Window Information.

Следует обратить внимание на три главных пункта меню:

File содержит подпункт New WinTool, с помощью вызова которого создаётся новое окно Window Design & Analysis Tool. Подпункт Full View Analysis открывает окно GUI Window Visualization Tool…, предназначенное для просмотра и экспорта окна в графический файл, как в растровый, так и в векторный. Пункт Export… (Рисунок 1.6) открывает диалоговое окно для экспорта данных фильтра в среду MATLAB, текстовый файл или матлаб-файл с возможностью выбора способа экспорта (коэффициентами или объектами) и возможностью изменить имя переменной, в которую экспортируется данные; для повторного экспорта в уже имеющуюся переменную необходимо поставить флажок Overwrite Variables, переписывать переменные.

Рис. 1.6. Диалоговое окно WinTool Export.

Меню View содержит подпункты, вызывая которые можно включать или отключать отображение графиков на первой панели, а также легенду графика частот. Важен подпункт Analysis Parameters…, позволяющий задать нужную степень аппроксимации функции частот (Number of points), развёртку (Range), величины по оси ординат (Response), единицы частоты (Frequency Units), вид графика (линейный/логарифмический, Frequency Scale).
Меню Tools (Рисунок 1.7), как и панель инструментов (Рисунок 1.8), находящаяся под главным меню, содержит инструменты для более подробного рассмотрения графиков. На графике можно поставить точки и наблюдать за изменениями функции (Рисунок 1.9).

Рис. 1.7. Меню Tools в окне Window Design & Analysis Tool.

Signal Processing Toolbox: Window Design and Analysis Tool - Menu and Panel

Рис 1.8. Главное меню и панель инструментов Window Design & Analysis Tool.

Рассмотрим функции меню Tools подробнее:

· Zoom In, увеличить, позволяет увеличить масштаб какого-либо из двух графиков обычным выделением нужной площади на графике (Рисунок 1.9), подобно выделению ярлыков, папок и файлов, удерживая левую кнопку мыши и двигая её, на Рабочем столе Windows. Такое же действие имеет кнопка панели инструментов:

Signal Processing Toolbox: Window Design and Analysis Tool - точки на графике

Рис. 1.9. Использование инструмента Zoom in из меню Tools. Про точки на графике и их использование будет сказано чуть ниже.

· Zoom x-axis, увеличить по оси абсцисс, растягивает график по оси абсцисс, сохраняя масштаб по оси ординат (Рисунок 1.10). Такое же действие имеет кнопка панели инструментов:

Signal Processing Toolbox: Window Design and Analysis Tool - Zoom x-axis
Рис. 1.10. Использование инструмента Zoom x-axis из меню Tools. В качестве примера взято окно Tukey. Его частотная характеристика увеличена по оси абсцисс (оси x) от 0 до 0.4 единиц, вместо стандартного масштаба от 0 до 1.
Используя этот инструмент, мышью надо создать прямую, параллельную оси абсцисс. Масштаб на оси x изменится в соответствии с положениями крайней левой и крайней правой точек этой прямой. Масштаб по оси y не изменяется.

· Zoom y-axis, увеличить по оси ординат, растягивает график по оси ординат, сохраняя масштаб по оси абсцисс. Действие инструмента подобно действию Zoom x-axis: создаётся прямая, параллельная оси ординат, причём верхняя и нижняя точка этой прямой указывают новую верхнюю точку и новую нижнюю точку всей оси соответственно. Аналог на панели инструментов:

· Full View, просмотреть графики полностью, уменьшает масштаб обоих графиков, если он был ранее увеличен и весгда растягивает временной график по оси ординат так, чтобы единица была максимальным значением оси ординат (оси Amplitude).

· Кнопка панели инструментов Restore Default View, восстановить отображение графиков по умолчанию, восстанавливает графики, если с ними были проведены какие-либо действия, связанные с изменением масштаба. Действие кнопки аналогично действию пункта меню Tools Full View за исключением растягивания графика временной характеристики. Кнопка имеет вид:

Окно Window Design & Analysis Tool даёт возможность с большой точностью наблюдать за координатами кривых. Например, точно узнаем координаты горбов частотной характеристики окна Тьюки (Tukey). Для этого нужно поместить указатель мыши на интересующую нас точку графика (можно делать это неточно) и щёлкнуть левую или правую кнопку мыши. Точки сохраняются при изменении масштабов, что очень удобно для отслеживания координат. Точки можно передвигать, причём одновременно отображаются числовые данные конкретной точки. Таким образом, можно обнаруживать максимумы (Рисунок 1.11) и минимумы горбов.

Signal Processing Toolbox: Window Design and Analysis Tool - нахождение максимумов горбов частотной характеристики окна Тьюки

Рис. 1.11. Нахождение максимумов горбов частотной характеристики окна Тьюки.

Для каждой точки, нажатием правой кнопки мыши, можно вызвать контекстное меню (Рисунок 1.12).

Signal Processing Toolbox: Window Design and Analysis Tool - контекстное меню для точки

Рис. 1.12. Контекстное меню для точки на частотном графике окна.

Рассмотрим содержимое контекстного меню (важные пункты выделены жирным):

Alignment, выравнивание, настраивает выравнивание данных о координатах относительно данной точки (справа сверху, слева сверху, справа снизу, слева снизу).
FontSize, размер шрифта, задаёт размер шрифта для данных о координатах.
Флаг Movable, возможность перемещения точки, показывает, можно ли двигать точку по графику или нельзя.
Interpolation, интерполяция. MATLAB строит графики по количеству точек, указанному в диалоговом окне Analysis Parameters (по умолчанию, 512), которое можно вызвать из меню View. Любой график, построенный на компьютере, содержит конечное количество этих точек, которое определяет точность построения графика. Для точки возможны два вида интерполяции: Nearest (Ближайшая точка) и Linear (Линейная интерполяция). При интерполяции Nearest, строго говоря, интерполяция в процессе работы программы не происходит – координаты конкретной точки берутся из уже сделанных программой отсчётов. Если рассмотреть график поближе, увидим, что в этом режиме невозможно расположить точку между отсчётами графика (вне точек излома), а сам график, на самом деле, представляет собой ломаную линию. Один из горбов графика, сильно увеличенный, показан на Рисунке 1.13.

Рис 1.13. Interpolation в режиме Nearest.
При интерполяции Linear возможно поставить точку на отрезок, соединяющий две соседние точки отсчёта, но значения абсциссы и ординаты точки будут неточными, так как интерполяция линейная.
Track Mode, режим передвижения, позволяет задавать привязку передвижения точки к передвижению указателя мыши.
Delete, Delete All удаляют точку или все точки на обоих графиках соответственно.

Полезно знать зависимости размерностей графиков по оси абсцисс:

Frequency Domain	Time Domain
Hz	sec
kHz	ms
MHz	µs
GHz	picosec

Если в диалоговом окне Analysis Parameters (меню View) в поле Frequency Units (единицы частоты) находится значение Normalized, то на оси x графика Time domain откладывается значение Samples. Это значит, что у окна нет конкретной длительности. Она может выражаться как в секундах, так и в миллисекундах, так и в микросекундах и т.д. Частота обратно пропорциональна времени, поэтому секундам на первом графике соответствуют герцы на втором, миллисекундам соответствуют килогерцы и т.д..

5. Window Visualization Tool

Окно Window Visualization Tool вызывается одной из команд:

>> wvtool(winname(n))
>> wvtool(winname₁(n),winname₂(n),...winname_m(n))
>> wvtool(w)

Где winname – имя окна в среде MATLAB, любая из оконных функций (см. параграф 5),n – количество отсчётов окна, то же, что и параметрL. w – вектор-столбец амплитуд (по нему строится график временной характеристики). Если в параметрах функции указано несколько оконных функций, то их графики будут выделены разным цветом. Рисунок 1.14 показывает результат выполнения команды построения треугольного окна из 64 отсчётов и окна Бартлетта из 64 отсчётов:

>> wvtool(triang(64),bartlett(64))

Signal Processing Toolbox: Window Visualization Tool

Рис. 1.14. Треугольное окно и окно Бартлетта в WVTool

Окно Window Visualization Tool создано для просмотра, лёгкого изменения и экспорта графиков окон. Меню Insert позволяет вставить на график прямоугольник (Rectangle), стрелку (Arrow), двустороннюю стрелку (DoubleArrow), линию (Line) и надпись (Textbox). Команда Export меню File позволяет сохранять изображения графиков в 15 форматах графических файлов, как в растровых, так и в векторных. Многие команды напоминают команды окна Window Design & Analysis Tools.

6. Виды окон и их характеристики (оконные функции)

Bartlett/Окно Бартлетта

w = bartlett(L)

Возвращает вектор-столбец отсчётов окна Bartlett в зависимости от параметра L. Параметр L – число отсчётов окна. L должен быть целым, большим нуля.

Отсчёты окна рассчитываются по формуле:

Причём L=N+1

Особенности:

Окно Bartlett очень похоже на треугольное окно. На концах окно всегда имеет нули, чего нет у треугольного окна. Центральные L-2 отсчёта окна Bartlett эквивалентны треугольным L-2 отсчётам (Рисунок 1.14, window #2).

Bartlett-Hanning/Окно Бартлетта-Хеннинга

w = barthannwin(L)

Особенности:

График подобен графику окна Bartlett. Как и для окон Bartlett, Hann и Hamming, частотная характеристика имеет асимптотически спадающие боковые лепестки (Рисунок 1.15). Это линейная комбинация взвешенных окон Бартлетта и Ханна c боковыми лепестками, меньше, чем у окон Бартлетта и Ханна вместе взятых. Ширина главного лепестка (mainlobe) не увеличилась, по сравнению с окнами Бартлетта и Ханна.

Отсчёты:

Где , причём .

Окна Бартлетта-Хеннинга, Бартлетта, Хемминга, Ханна.

Рис. 1.15. #1 - Bartlett-Hanning, #2 - Bartlett, #3 – Hamming, #4 – Hann.

Blackman/Окно Блэкмана

w = blackman(L)
Отсчёты:

L=N+1

Blackman-Harris/Окно Блэкмана-Харриса

w = blackmanharris(L)

Особенности:

Окно оптимизировано по площади боковых лепестков. В формуле вычисления коэффициентов - 4 слагаемых.

Отсчёты:

Bohman/Окно Бохмана

w = bohmanwin(L)

Отсчёты:

x – вектор L, но линейно разделённый (необходимо использовать linspace).

Термин «линейно разделённый» означает, что разница предыдущего и текущего элементов вектора равна разнице текущего и последующего элементов.

Особенности:

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

Chebyshev/Окно Чебышева

w = chebwin(L,r)

r – параметр, который показывает, на сколько дБ главный лепесток выше боковых лепестков. По умолчанию этот параметр равен 100 дБ.

Особенности:

Присутствие импульсов в конечных точках отсчёта временной характеристики. Это вызвано постоянным уровнем боковых лепестков в частотной характеристике. Если величина боковых лепестков большая, эффект в крайних точках временной характеристики может быть существенным.

Рассмотрим выполнение команды:

>> wvtool(chebwin(64,20),chebwin(64,40),chebwin(64,100),chebwin(64,200))

Окна Чебышева

Рис. 1.16. Функиция chebwin. #1 – r=20, #2 – r=40, #3 – r=100, #4 – r = 200.

Уровень боковых лепестков эквивалентен задержке в полосе подавления. Исходя из графика 1.16, чем ниже уровень боковых лепестков, тем шире полоса пропускания, тем более полого и неидеально спадает частотная характеристика фильтра.

Увеличение ширины главного лепестка АЧХ окна, при той же длительности импульсной характеристики приводит к «размазыванию» АЧХ фильтра - крутизна АЧХ в области частоты среза уменьшается (расширяется переходная полоса). Чтобы сохранить крутизну АЧХ, необходимо увеличивать длительность импульсной характеристики фильтра.

Flat Top/Окно с плоской полосой пропускания

w = flattopwin(L)

Особенности:

Окно Flat Top имеет очень маленькие пульсации в полосе пропускания (менее 0.01 дБ). Используется для калибровки. Полоса пропускания приблизительно в 2.5 раза шире, чем у окна Ханна.

Отсчёты:

, везде кроме n.

L=N+1

a₀= 0.21557895 a₁= 0.41663158 a₂= 0.277263158 a₃= 0.083578947 a₄= 0.006947368.

Gaussian/Окно Гаусса

w = gausswin(N)
w = gausswin(N,Alpha)

N - число отсчётов. N=L.

Коэффициенты:

Особенности:

Ширина окна обратно пропорциональна параметру альфа (Рисунок 1.17). Больший параметр альфа создаёт более узкое окно. По умолчанию альфа принимается равным 2.5.

Окна Гаусса

Рис 1.17. Окно Gaussian. #1 – α=1, #2 – α=2, #3 – α =2.5, #4 – α =5.

Hamming/Окно Хемминга

Графики изображены на Рисунке 1.15

w = hamming(L)

Отсчёты:

L=N+1.

Hann/Окно Ханна

Графики изображены на Рисунке 1.15

w = hann(L)

Отсчёты:

L=N+1

Kaiser/Окно Кайзера

w = kaiser(L,beta)

Параметр бета указывает уровень ослабления боковых лепестков.

Для создания окна Кайзера с ослаблением лепестковαдБ, нужно использовать следующие значения β:

>> wvtool(kaiser(200,0),kaiser(200,3.39532),kaiser(200,10.06126)) %Рис. 1.18

Окна Кайзера

Рис. 1.18. Окно Kaiser. #1 – α<21, β=0; #2 – α=40,β=3.39532; #3 – α=100, β=10.06126.

Как и для других окон, при изменении боковых лепестков в желаемую сторону (на графике – ниже, по физическому смыслу – меньше пропускания в полосе подавления), расширяется полоса пропускания и становится более пологим переход от полосы пропускания к полосе подавления. Более подробно график частот можно рассмотреть на Рисунке 1.19.

Окна Кайзера

Рис. 1.19. Частотная характеристика окна Kaiser в зависимости от коэффициентов.

Nuttall

w = nuttallwin(L)

Боковые лепестки оптимизированы по мощности. Окно появилось на основе окна Блэкмана-Харриса.

>> L = 64;

w = blackmanharris(L);

y = nuttallwin(L);

wvtool(w,y) % Рис. 1.20

Окна Натталла и Блэкмана-Харриса

Рис. 1.20. Окна Nuttall (#2) и Blackman-Harris(#1).

Окно Nuttall имеет боковые лепестки на одном уровне.

Отсчёты:

a₀ = 0.3635819, a₁ = 0.4891775, a₂ = 0.1365995, a₃ = 0.0106411

Parzen/Окно Парзена

w = parzenwin(L)

Окна Parzen – кусочные кубические аппроксимации окна Gaussian. Боковые лепестки окна Parzen убывают как 1/w⁴.

Отсчёты:

На интервале:

>> wvtool(parzenwin(64),gausswin(64)) % Рис. 1.21.

Окно Парзена и Гаусса

Рис. 1.21. Окна Parzen (синий, оба графика снизу) и Gaussian (зелёный).

Rectangular/Прямоугольное окно

w = rectwin(L)

Отсчёты рассчитываются с помощью функции:

w = ones(L,1);
Y = ones(m,n) – формирует массив единиц Y размерности m на n.

Triangular/Треугольное окно

triang(L)Окно и его частотная характеристика изображены на Рисунке 1.14, window #1.Отсчёты:Для нечётных:

Для чётных:

Tukey/Окно Тьюки

w = tukeywin(L,r)Особенности:Окно Тьюки с обеих сторон представляет собой функцию косинуса.Коэффициент r равен отношению постоянной составляющей окна к составляющей в виде функции косинуса. r не может быть меньше нуля или больше единицы.
Отсчёты:

Причём α=r. Если указать α≤0, создаётся прямоугольное окно, если указать α≥1, создаётся окно Ханна. Окно Тьюки

Рис. 1.22. Окно Tukey при различных параметрах r. #1:r=0, #2:r=0.5, #3:r=0.8.

Некоторые выводы

Чем ниже уровень боковых лепестков, тем шире полоса пропускания и тем больше переходная характеристика отличается от идеальной. Известно, что увеличение ширины главного лепестка окна при той же длительности импульсной характеристики приводит к «размыванию» АЧХ (Frequency domain, частотной характеристики) фильтра – крутизна АЧХ в области частоты среза уменьшается (расширяется переходная полоса). Чтобы сохранить крутизну АЧХ, необходимо увеличивать длительность импульсной характеристики.
Можно выделить группы окон с некоторыми похожими свойствами:
1. Nuttall и Blackman-Harris: одинаковые полосы пропускания, разные полосы подавления. У Nuttall полоса подавления сразу после главного лепестка ниже, чем в других точках частотной характеристики, а у окна Blackman-Harris наоборот: после главного лепестка боковые лепестки выше, чем в других точках частотной характеристики. (Рисунок 1.20.)
2. Окна Kaiser, Gaussian, Chebyshev при различных параметрах. В зависимости от своих параметров, полоса перехода различна, но в пределах одного окна поведение лепестков в зависимости от частоты является одинаковым. Так, в окнах Chebyshev (Рисунок 1.16) и Gaussian (Рисунок 1.17) уровень лепестков одинаков, а в окне Kaiser (Рисунок 1.18) наблюдается небольшое падение лепестков и их асимптотическое стремление.
3. У окон Bartlett-Hanning, Bartlett, Hamming, Hann (Рисунок 1.15) одинакова полоса пропускания, но сильно различаются боковые лепестки.
4. У окон Bartlett и Triangular (Рисунок 1.14) различия на первых лепестках малы (различие в 10 дБ – только на седьмом боковом лепестке), но постепенно возрастают.

7. Перемножение характеристики идеального фильтра на характеристику окна

Используя команду

b = 0.4*sinc(0.4*(-25:25));
создадим идеальный КИХ-фильтр нижних частот с 51 отсчётом (Рисунок 1.23). Частота среза фильтра – 0.4π.
Функция sinc определяется как:

Функция sinc является обратным преобразованием Фурье от прямоугольного импульса ширины 2π и высоты 1:

Рис 1.23. Идеальный КИХ-фильтр.
Умножим вектор отсчётов идеального фильтра b на окно Hamming (Рисунок 1.24):
b = b.*hamming(51)';% .* - операция умножения матриц,% ' – транспонирование.
АЧХ после перемножения sinc на окно

Рис 1.24. Перемножение на окно. Результирующая АЧХ.
Так происходит перемножение на окно.

8. Функции fir1 и fir2

Функция fir1 имеет следующие параметры (представлены основные возможности):
b = fir1(n,Wn)

b = fir1(n,Wn,'ftype')

b = fir1(n,Wn,window)

b = fir1(n,Wn,'ftype',window)
n – порядок фильтра.
Wn – частота среза фильтра, может быть как одной переменной, так и вектором. Если Wn – фектор, и Wn = [w1 w2 w3 w4 w5 ... wn], то частоты пропускания определяются следующим образом: 0 < ω< w1, w1 < ω< w2, ..., wn < ω< 1, создаётся многополосный фильтр.
Параметр 'ftype' определяет тип фильтра: 'high' – для фильтра верхних частот, 'stop' – для создания режекторного фильтра, причём Wn должна быть двухэлементным вектором, который указывает границы полосы задержки, 'DC-1' – для создания в многополосном фильтре первой полосы в качестве полосы пропускания, 'DC-0' – для создания первой полосы в многополосном фильтре в качестве полосы подавления.
Параметр window определяет окно, причём длина окна должна быть равна n+1. Если параметр не указан, проектирование проходит с использованием окна Hamming.
Групповое время задержки фильтра, спроектированного с помощью fir1, равно n/2.
Функция fir2 имеет следующие параметры (представлены основные возможности):
b = fir2(n,f,m)

b = fir2(n,f,m,window)
f - вектор значений частоты.
m – вектор значений амплитуды при данных частотах f, должно соблюдаться равенство: длина вектора f = длине вектора m.

9. Объекты sigwin

При проектировании фильтра с помощью окон могут оказаться полезными объекты sigwin с полями данных Name (имя) и Length (количество отсчётов).

Операция создания объекта типа окно:

w=sigwin.имя_окна(параметры окна)

При отсутствии параметров будет создано окно, содержащее 64 отсчёта.

Для объекта sigwin используются следующие функции:

get(w) % вывод всех значений полей окна
get(w,'parameter') % вывод значения указанного параметра окна
set(w,'parameter1',value,'parameter2',value,...) % изменение значений полей
vector=generate(w) % создание вектора отсчётов vector окна w
info(w) % вывод информации об окне
w.winwrite('имя_файла') % вывод информации об окне w в файл
% имя_файла.wf в папке Мои Документы/MATLAB/
Пример использования:
>> w=sigwin.bartlett; % создан объект окна с 64 отсчётами
>> get(w,'Length') % вывод значения параметра Length
ans = 64
>> set(w,'Length',128) % изменение значения параметра Length
>> info(w) % вывод всей информации об окне
ans =Bartlett Window
---------------
Length : 128
>> v=generate(w); % создание вектора отсчётов окна
>> w.winwrite('file'); % вывод отсчётов и информации об окне в файл file.wf

Критерии точности фильтрации сигнала КИХ-фильтром

1. Часто используемые критерии точности

Существует множество критериев точности фильтрации:· Разрядность сигнальной памяти· Разрядность арифметических действий· Разрядность АЦП· Устойчивость при определённом сигнале· Дисперсия, максимальная разность и другие математические методы

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

2. Специфическая оценка окон

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

Equivalent noise bandwidth / Эквивалентная шумовая полоса

Из Рисунка 2.1 видно, что оценка амплитуды гармонической компоненты на заданной частоте оказывается смещенной из-за наличия широкополосного шума, попадающего в полосу пропускания окна. К энергии полезного сигнала добавляется энергия шума, ослабленного окном. (Шумом является всё, что не является нужной гармоникой, в том числе «лепестки» от этой гармоники.) В этом смысле окно ведет себя как фильтр, мощность сигнала на выходе которого пропорциональна мощности гармоник входного сигнала в полосе его пропускания. Для обнаружения гармонического сигнала необходимо минимизировать накопленный шум. Этого можно достичь с помощью узкополосного окна. Идеальное окно в частотной области представляет собой дельта-функцию, однако это окно нереализуемо. Удобной мерой ширины полосы пропускания окна является его эквивалентная шумовая полоса (ЭШП, ENBW в английской литературе). ЭШП окна – это ширина полосы пропускания прямоугольного окна с тем же максимальным усилением по мощности, который накапливает ту же мощность шума, что и данное окно (Рисунок 2.2). Другими словами, эквивалентная шумовая полоса – это число, большее единицы и показывающее, какой должна быть ширина полосы пропускания прямоугольного окна, если оно будет накапливать ту же мощность шума, что и данное окно, причём это прямоугольное окно во временной области должно быть такой же высоты в частотной области, как и максимальная высота оцениваемого окна.

Накопленная окном мощность шума (Noise Power) определяется выражением:

(Формула 2.1)

N₀ – мощность шума на единицу полосы пропускания.

Согласно теореме Парсеваля, эту величину можно выразить следующим образом:

(Формула 2.2)

Рис. 2.1. Иллюстрация свёртки. Окно представлено в виде спектрального фильтра.

Эквивалентная шумовая полоса окна

Рис. 2.2. Эквивалентная шумовая полоса окна.

Максимальное усиление по сигналу (Peak Signal Gain) и по мощности (Peak Power Gain) соответствует частоте ω=0, оно называется усилением по сигналу / по мощности на нулевой частоте и определяется выражениями:

(Формула 2.3)

ЭШП (ENBW в английской литературе) окна, нормированная на величину N₀/T – мощность шума на единичный временной интервал, может быть записана в виде:

(Формула 2.4)

Таким образом, эквивалентная шумовая полоса окна прямо пропорциональна (с коэффициентом пропорциональности из-за нормирования N₀/T) мощности шума, которую накопило окно, и обратно пропорциональна максимальному усилению окна по мощности.

Processing Gain / Усиление преобразования

С ЭШП окна тесно связаны понятия усиления преобразования (УП, PG, Processing Gain в английской литературе) и потерь преобразования (ПП, PL, Processing Lost в английской литературе) при вычислении ДПФ с помощью окон. ДПФ можно рассматривать как результат пропускания сигнала через набор согласованных фильтров, каждый из которых настроен на одну из гармоник комплексной синусоидальной последовательности базисного множества, состоящего из функций косинуса и синуса. С этой точки зрения мы и будем анализировать усиление преобразования (называемое также когерентным усилением) фильтра и потери преобразования, вызванные тем, что окно сглаживает, т. е. сводит к нулю, величины отсчетов, расположенных вблизи его границ. Пусть входная последовательность отсчетов задана выражением:

(Формула 2.5)

q(nT) – последовательность отсчетов белого шума с дисперсией .

Составляющая сигнала в спектре, вычисленном с помощью окна, будет равна:

(Формула 2.6)

Из (2.6) видно, что в отсутствие шума спектральная составляющая пропорциональна входной амплитуде А. Таково же будет и математическое ожидание этой составляющей при наличии шума. Коэффициент пропорциональности равен сумме всех отсчетов дискретного окна, а эта сумма есть не что иное, как усиление окна для постоянного сигнала. Для прямоугольного окна этот коэффициент равен N - числу отсчетов в окне. Усиление любого другого окна меньше, так как весовая функция вблизи границ окна плавно спадает к нулю. Связанное с этим уменьшение коэффициента пропорциональности полезно знать, так как оно характеризует ошибку (смещение) оценок амплитуд спектральных составляющих. В литературе вместо усиления преобразования иногда используется другой параметр – когерентное усиление по мощности, т. е. квадрат когерентного усиления сигнала.

Некогерентная составляющая взвешенного, т. е. выполненного с помощью окна преобразования, вычисляется по формуле:

(Формула 2.7)

Некогерентная мощность определяется выражением:

(Формула 2.8)

E {} – оператор математического ожидания.

Заметим, что некогерентное усиление по мощности равно сумме квадратов отсчетов весовой функции, а когерентное – квадрату суммы этих отсчетов.

И наконец, вычислим усиление преобразования (УП, PG в английской литературе), которое определяется как частное от деления отношений сигнал/шум на выходе и на входе:

(Формула 2.9)

Заметим, что усиление преобразования – это величина, обратная нормированной ЭШП окна. Таким образом, увеличение ЭШП окна ведет к уменьшению УП. Это вполне понятно, так как, чем шире полоса пропускания, тем больше мощность шума, прошедшего через окно и вносящего вклад в спектральную оценку.

Overlap Correlation / Корреляция перекрываемых участков

Когда быстрое преобразование Фурье (БПФ) используется для обработки длинной последовательности, эту последовательность предварительно делят на несколько последовательностей по N отсчетов каждая, при этом N выбирается так, чтобы обеспечить необходимое спектральное разрешение. Спектральное разрешение БПФ определяется Формулой 2.10, где ∆f – спектральное разрешение, f_s – частота дискретизации, выбранная согласно критерию Найквиста, и β – коэффициент, характеризующий увеличение ширины полосы для выбранного окна. Отметим, что f_s/N – это наилучшее разрешение, достижимое при БПФ. Коэффициент β обычно выбирается равным ЭШП окна в бинах (бин – единичный временной интервал).

(Формула 2.10)

Если окно и БПФ «воздействуют» на неперекрывающиеся участки последовательности (Рисунок 2.3), то значительная часть данных попросту игнорируется, поскольку вблизи границ окна значения его отсчетов близки к нулю. Так, например, если преобразование используется для обнаружения коротких узкополосных сигналов, то при анализе неперекрывающихся участков появление сигнала может оказаться просто незамеченным. Для этого достаточно, чтобы сигнал появился вблизи границы любого из интервалов. Чтобы избежать таких потерь данных, преобразованию обычно подвергают перекрывающиеся участки последовательности. Степень перекрытия в большинстве случаев выбирается равной 50 или 75%. Разбиение сигнала на перекрывающиеся участки, конечно, увеличивает общий объем вычислений, однако достигаемые с его помощью результаты вполне это оправдывают.

Перекрытие

Рис. 2.3. Перекрытие последовательностей.

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

(Формула 2.11)

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

(Формула 2.12)

Зададимся теперь вопросом, насколько уменьшится дисперсия при усреднении коррелированных измерений, как это имеет место при усреднении преобразований Фурье перекрывающихся участков? Ответ на это дал Уэлш. Результат без доказательств для случая 50% и 75% перекрытия:

(Формула 2.13)

Отрицательные члены в (2.13) описывают краевые эффекты усреднения; их можно не учитывать при K>10. Для хороших окон член c²(0,25)<<1, и его также можно опустить с пренебрежимо малой ошибкой. Заметим, что для хороших окон преобразование перекрывающихся на 50% участков сигнала практически независимы.

Scalloping Loss / Паразитная амплитудная модуляция спектра

Важным фактором, влияющим на обнаружимость слабых сигналов, является паразитная амплитудная модуляция спектра (scalloping loss),или эффект «частокола» (picket-fence effect). Ранее мы рассматривали выполняемое с помощью окна ДПФ как результат пропускания сигнала через набор согласованных фильтров и анализировали обусловленные специфическими свойствами окна усиление и потери для тонов, совпадающих с базисными векторами. Базисные векторы – это тоны, кратные частоте f_s/N, где f_s – частота отсчётов. Эти частоты не что иное, как точки отсчетов спектра, их обычно называют выходными точками, частотами гармоник или бинами ДПФ. Зададимся теперь вопросом, каковы будут дополнительные потери при обработке сигнала, частота которого лежит посредине между частотами соседних бинов (т. е. сигнала с частотой (k+l/2)f_s/N)?

Обратившись к формуле 2.6 и заменив в ней ω_k на ω_k_+1/2 получаем, что усиление окна для частоты, сдвинутой на 0.5 бина, равно:

(формула 2.14)

По определению, потери из-за паразитной амплитудной модуляции (AM, Scallop Loss, Scalloping Loss, SL в английской литературе) спектра равны отношению когерентного усиления тона, расположенного посредине между двумя бинами (бин – единичный временной интервал) ДПФ, к когерентному усилению тона, совпадающего с одним из бинов ДПФ, то есть:

(Формула 2.15)

Потери из-за паразитной амплитудной модуляции равны максимальным потерям при самой неблагоприятной для ДПФ частоте сигнала.

Worst Case Processing Lost / Максимальные потери преобразования

Теперь сделаем одно интересное замечание. Определим максимальные потери преобразования (ПП) как сумму максимальных потерь из-за паразитной амплитудной модуляции спектра для данного окна (в дБ) и потерь преобразования, обусловленных формой этого окна. Введенный параметр характеризует уменьшение выходного соотношения сигнал/шум в результате воздействия окна при наихудшем расположении частоты сигнала. Его величина, конечно, влияет на минимальную интенсивность тона, при которой он еще может быть обнаружен в широкополосном шуме. Интересно заметить, что уровень максимальных потерь всегда лежит между 3,0 и 4,3 дБ. Окна, для которых максимальные ПП превышают 3,8 дБ, совершенно неудовлетворительны, и их не следует применять. Из сопоставления данных о потерях, приведенных в Таблице 1 и на Рисунке 2.6, видно, что почти все окна (за исключением прямоугольного) одинаково пригодны для обнаружения чистых тонов в широкополосном шуме. Разница в потерях у различных окон не превышает 1,0 дБ, а для хороших окон – 0,7 дБ. Однако обнаружение тона в присутствии других близких тонов представляет собой совершенно иную задачу. Именно здесь тип применяемого окна может иметь решающее значение.

Ещё раз о просачивании спектральных составляющих

Возвращаясь к свёртке и Рисунку 2.1, заметим, что на точность измерения амплитуды спектральной составляющей влияет не только спектр широкополосного шума, но и узкополосные помехи, если они попадают в полосу пропускания окна. Действительно, некоторая спектральная компонента, скажем, с частотой ω=ω₀ будет вносить вклад в спектральную компоненту с частотой ω= ω_a, т. е. будет наблюдаться на этой частоте. Этот вклад будет определяться усилением окна с центром в ω₀ на частоте ω_a. Это и есть эффект, называемый просачиванием спектральных составляющих. Он показан на Рисунке 2.4 для преобразования конечного тона частотой ω₀.

Просачивание приводит к смещению оценок амплитуд и положений гармонических составляющих сигнала. Даже для единственной вещественной гармоники (не совпадающей с частотой гармоники ДПФ) просачивание от ядра на оси отрицательных частот влияет на ядро на оси положительных частот. Это влияние наиболее сильно и неприятно при обнаружении слабых сигналов в присутствии сильных помех близкой частоты. Для уменьшения неприятных последствий из-за спектрального просачивания амплитуда боковых лепестков вдали от главного центрального лепестка частотной характеристики окна должна быть малой, а переход от центрального лепестка к низкоамплитудным боковым лепесткам – очень быстрым. Одним из параметров, указывающих, насколько хорошо окно подавляет просачивание, является максимальный уровень боковых лепестков (по отношению к главному лепестку), другой параметр – это асимптотическая скорость спада боковых лепестков.

Просачивание спектра

Рис 2.4. Просачивание спектральных составляющих.

Minimum Resolution Bandwidth / Минимальная разрешаемая полоса частот

Спектральное разрешение

Рис. 2.5. Спектральное разрешение двух близкорасположенных ядер.

Рисунок 2.5 подсказывает еще один критерий, который должен использоваться при выборе оптимальных окон. Поскольку окно придает спектральной линии некоторую эффективную ширину, интересно знать, при каком минимальном расстоянии между двумя спектральными линиями равной интенсивности главные лепестки этих линий еще могут быть разрешены независимо от положения линий относительно бинов ДПФ. Классический критерий такого разрешения – ширина окна между точками, в которых мощность главного лепестка спадает наполовину (ширина окна по уровню 3,0 дБ). Этот критерий отражает тот факт, что два главных лепестка равной интенсивности, отстоящие друг от друга по частоте менее чем на ширину окна по уровню 3,0 дБ, будут иметь один общий спектральный пик и не будут разрешаться как две отдельные линии. Однако трудность использования этого критерия в том, что он несовместим с когерентным суммированием, применяемым в ДПФ. Выходные точки ДПФ получаются путем когерентного сложения спектральных компонент, взвешенных окном с центром на данной частоте.

Если в когерентное суммирование вносят вклад два ядра, их сумма в точке пересечения (номинально посредине между ними) должна быть меньше, чем индивидуальные пики, если эти пики должны быть разрешены. Таким образом, в точках пересечения ядер усиление от каждого ядра не должно превышать 0,5, т. е. расстояние между пиками должно превышать ширину окна по уровню 6,0 дБ. Ширина различных окон по уровню 6,0 дБ указана в Таблице 1. Из таблицы видно, что эта ширина меняется от 1,2 до 2,6 бин, где бин – разрешение основной частоты ω_s/N. Что касается ширины окна по уровню 3,0 дБ, то, как будет показано в следующем абзаце, она также является полезной характеристикой окна. Следует, однако, помнить, что разрешение ДПФ определяется шириной используемого окна по уровню 6,0 дБ.

Некоторые итоги и классификация

Из Таблицы 1 видно, что шумовая полоса (ЭШП, ENBW) окна всегда шире его полосы по уровню 3,0 дБ (3.0-dB BW). Разность этих двух параметров, отнесенная к ширине окна по уровню 3,0 дБ, является довольно чувствительным показателем качества работы окна. Для всех хороших окон, перечисленных в Таблице 1, этот показатель лежит в пределах 4,0–5,5%. Окна, для которых он лежит вне этого диапазона, имеют либо широкий главный лепесток, либо высокий уровень боковых лепестков. Для таких окон характерны либо высокие потери преобразования, либо малая эффективность обнаружения двух близких тонов. Окна, для которых этот показатель лежит в пределах 4,0–5,5%, попадают в левый нижний угол диаграммы на Рисунке 2.6, характеризующей качество их работы.

В Таблице 1 перечислены общепринятые параметры, характеризующие типы окон. Однако ориентироваться в этой массе цифр непросто. Заметим, что самым важным из табличных параметров, вероятно, является уровень боковых лепестков (чем он ниже, тем меньше смещение спектральных оценок). Другой важнейший параметр – максимальные потери преобразования (чем они ниже, тем выше обнаружимость слабых сигналов). Рисунок 2.6 показывает относительное положение различных окон как функцию этих двух параметров. Точки, соответствующие лучшим окнам, лежат в левом нижнем углу диаграммы. Такие окна имеют низкий уровень боковых лепестков и низкие максимальные потери преобразования.

Сравнительная характеристика окон как функций от высочайшего бокового лепестка и максимальных потерь преобразования

Рис 2.6. Классификация окон по максимальным потерям и наибольшему боковому лепестку. Лучшие окна лежат в левой нижней части.

Влияние боковых лепестков

Из вышесказанного отметим вкратце важные параметры для окон:

Эквивалентная шумовая полоса (ЭШП, ENBW) зависит от мощности шума и усиления окна, которые, в свою очередь, вычисляются посредством использования временной области и операций с её отсчётами.
Усиление преобразования (УП, PG), частное от отношения сигнал/шум на выходе и на входе, зависит, также как и ЭШП, от величин отсчётов на временной области.
Корреляция перекрываемых участков – давно устоявшаяся теория перекрытия окон, далее рассматривать её мы не будем.
Паразитная амплитудная модуляция спектра (Scallop Loss), т.е. дополнительные потери при обработке сигнала, лежащего по частоте между частотами отсчётов ДПФ, также зависит от величин отсчётов на временной области.
Максимальные потери преобразования (Worst Case Process Lost) зависят от вида окна во временной области.
Просачивание спектральных составляющих зависит от лепестков спектра: максимального и уровня спада.
Минимальная разрешаемая полоса частот является характеристикой, не зависящей от боковых лепестков, она зависит от вида окна и далее рассматриваться не будет.

Временная область и частотная область связаны преобразованием Фурье и зависят друг от друга, поэтому разность ЭШП и полосы по уровню в 3 дБ, отнесённая к ширине окна по уровню в 3 дБ (важная характеристика) качественно может оцениваться по уровню главного лепестка и уровню боковых лепестков. Как отмечалось, это отношение влияет на обнаружимость двух близких тонов и на потери преобразования.

Итак, боковые лепестки (уровень главного лепестка, скорость спада) влияют на:

Просачивание спектральных составляющих
Обнаруживаемость двух близких тонов
Потери преобразования

Кроме того, Харрис отмечает, что «низкий уровень боковых лепестков соответствует лучшим окнам».

Оценим это влияние на практике.

Фильтрация разных видов сигналов тремя фильтрами: выбранным для взвешивания, спроектированным с помощью прямоугольного и спроектированным с помощью хорошего окна в MATLAB. Спектральные характеристики. Выводы.

1. Выбор фильтра для взвешивания

Соображения по выбору фильтра для последующего взвешивания

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

(Формула 3.1)

В MATLAB эту функцию можно описать так:

imp = [1; zeros(49,1)]; % Вектор, в котором одна единица и % 49 нулей, начиная со второй позиции, причём% количество нулей здесь очень важно:% длина вектора совпадает% с длиной вычисляемой импульсной характеристикиТогда импульсная характеристика фильтра:h = filter(b,a,imp); % b, a – коэффициенты числителя и знаменателя% функция filter фильтрует сигнал imp% через фильтр с коэффициентами b и a% если используется КИХ-фильтр, то a=1

Проектирование с помощью окон (windowing в английском) – метод, используемый при создании фильтра с конечной импульсной характеристикой. Кроме использования весовой функции, а окно – это и есть весовая функция, можно использовать ещё 4 метода проектирования КИХ-фильтра из пакета MATLAB. Некоторые из них основаны на понятии об ошибках и их минимизации, а некоторые даже могут синтезировать фильтр с нелинейной фазовой характеристикой (функция cfirpm)! Проектирование КИХ фильтра с помощью весовой функции в пакете MATLAB чаще всего осуществляется с помощью идеального фильтра с конечной импульсной характеристикой (функция sinc) и самого окна (как было ранее сказано, их импульсные характеристики перемножаются). Может и можно для этих целей брать любой неидеальный КИХ-фильтр, но в этом нет смысла, поскольку этот фильтр тоже спроектирован каким-либо другим способом.

Проектировать КИХ-фильтр с помощью взвешивания БИХ-фильтра сложнее, если такое, вообще говоря, возможно. Причина кроется в импульсной характеристике БИХ-фильтра (Рисунок 3.1). Импульсная характеристика КИХ-фильтра симметричная, как и отсчёты временной области любого из окон, совпадает с набором коэффициентов числителя передаточной функции b и показывает веса для элементов входного временного ряда (Рисунок 3.2). Импульсная характеристика БИХ-фильтра несимметричная, а передаточная функция имеет коэффициенты как в числителе, так и в знаменателе (Рисунок 3.1).

Импульсная характеристика БИХ-фильтра

Рис. 3.1. Импульсная характеристика БИХ-фильтра (слева).

Импульсная характеристика КИХ-фильтра

Рис. 3.2. Импульсная характеристика КИХ-фильтра (слева).

Если перемножить импульсную характеристику фильтра и импульсную характеристику окна, то будут ослаблены участки импульсной характеристики фильтра на краях. Стоит отметить, что участки импульсной характеристики фильтра посередине не получают усиления, так как максимальный отсчёт любого окна не превышает единицы. Принимая это во внимание, умножим импульсную характеристику БИХ-фильтра на окно (Рисунок 3.3) и получим соответствующую импульсную характеристику.

Умножение импульсной характеристики БИХ-фильтра на окно

Рис. 3.3. Результат умножения импульсной характеристики БИХ-фильтра на окно.

Естественно, хорошей АЧХ от такой операции не получилось. Напрямую импульсную характеристику БИХ-фильтра нельзя умножать на коэффициенты окна.

Стоит сказать, что иногда удобно применять GUI wvtool для отображения импульсной характеристики и АЧХ обычных фильтров (Рисунки 3.1, 3.2, 3.3). Если это КИХ-фильтр, то параметром wvtool служит вектор коэффициентов числителя, а если это БИХ-фильтр, то параметром wvtool может служить его импульсная характеристика (но ни в коем случае не коэффициенты!).

Выбор фильтра для последующего взвешивания

Для последующего взвешивания выберем функцию sinc, о которой уже рассказывалось в первой части работы. Эту функцию приводят в качестве примера для последующего взвешивания в Help MATLAB. Её используют чаще всего при проектировании с помощью окон.

b = 0.4*sinc(0.4*(-25:25));

Импульсная характеристика и АЧХ функции показаны на Рисунке 3.4 и Рисунке 3.5 соответственно.

Импульсная характеристика идеального фильтра, полученного функцией MATLAB sinc

Рис. 3.4. Импульсная характеристика 0.4*sinc(0.4*(-25:25))

АЧХ идеального фильтра, полученного функцией MATLAB sinc

Рис. 3.5. АЧХ 0.4*sinc(0.4*(-25:25))

Первый коэффициент (0.4) влияет на высшую точку импульсной характеристики, второй коэффициент (параметр функции sinc) влияет на частоту среза фильтра, третий и четвёртый коэффициенты (-25:25) влияют на количество отсчётов импульсной характеристики.

2. Взвешивание фильтра

Взвесим фильтр sinc с помощью плохого прямоугольного и хорошего, входящего в левую нижнюю четверть Рисунка 2.6, окна Блэкмана:

brect = b.*rectwin(51)'; % синий цвет на рисунке
bblack = b.*blackman(51)'; % зелёный цвет на рисунке

Результат показан на Рисунках 3.6 и 3.7.

Рис. 3.6. Спроектированные КИХ-фильтры. Амплитуда указана в децибелах.

Импульсная характеристика и АЧХ взвешенных фильтров

Рис. 3.7. Спроектированные КИХ-фильтры. Амплитуда указана на линейной шкале.

3. Выбор сигнала и пропускание его через фильтры.

Фильтрация синусоидального сигнала

Рис. 3.8. Фильтрация синусоидального сигнала.

b = 0.4*sinc(0.4*(-25:25)); % вычисление коэффициентов

% идеального фильтра

brect = b.*rectwin(51)'; % вычисление коэффициентов фильтра

% с помощью прясмоугольного окна

bblack = b.*blackman(51)'; % вычисление коэффициентов фильтра

% с помощью окна Блэкмана

t = 0:50; % задание времени, 51 отсчёт

x=cos(50*t); % вычисление косинусоиды из 51 отсчёта

y = filter(b,1,x); % фильтрация идеальным фильтром

yrect=filter(brect,1,x);

yblack=filter(bblack,1,x); % фильтрация неидеальным фильтром

xcoord=0:50; % задание координаты по оси абсцисс для графика

plot(xcoord,x,xcoord,y,xcoord,yrect,xcoord,yblack)

% data1 – сам сигнал, data2, data3 (слились) – фильтрация идеальным и

% фильтром, построенным с помощью прямоугольного окна (красный)

% data4 – фильтрация сигнала фильтром, спроектированным с помощью

% окна Блэкмана (бирюзовый)

Рисунок 3.8 и 3.9 показывает, что у фильтра, спроектированного с помощью окна Блэкмана гораздо лучше переходные характеристики. В фильтрах в начальный момент времени идёт переходной процесс, предполагается, что сигнал возник в нулевой момент времени, до этого его уровень был равен нулю. Длительность переходного процесса КИХ-фильтра равна половине его импульсной характеристики. После завершения переходного процесса видно, что фильтр задерживает входной сигнал. Максимальные значения колебаний переходного процесса равны, но видно, что фильтр, спроектированный с помощью прямоугольного окна, ведёт себя хуже в начальные моменты времени. Максимальная разница между отфильтрованными сигналами в эти моменты – 2%.

Переходные характеристики (первая половина отсчётов)

Рис. 3.9. Переходные характеристики разных фильтров. Красным цветом показана характеристика фильтра, спроектированного с помощью окна Блэкмана, зелёным – с помощью прямоугольного окна.

Также можно подтвердить выводы о том, что прямоугольное окно равносильно отсутствию окна – графики y и yrect совпадают.

При оценке дисперсии (в MATLAB эта оценка осуществляется с помощью функции var(x)) может оказаться, что дисперсия сигнала на выходе фильтра, спроектированного с помощью окна Блэкмана больше, чем дисперсия сигнала на выходе фильтра, спроектированного с помощью прямоугольного окна. Но это не должно служить показателем того, что окно Блэкмана хуже – на самом деле оно лучше передаёт амплитуды входного сигнала, что видно из графика (Рисунок 3.10). Ошибка фильтра, спроектированного с помощью прямоугольного окна – 3 сотых, ошибка хорошего фильтра при данном синусоидальном сигнале – 1 сотая (разница в 1% по отношению к удвоенной амплитуде входного сигнала в 2 единицы).

Рис. 3.10. Входной сигнал (слева, синий) и его передача различными фильтрами (справа, красный и бирюзовый). Цвета совпадают с цветами на Рисунке 3.8 (см. комментарии в коде).

Фильтрация экспоненциального сигнала

Рис. 3.11. Фильтрация экспоненциального сигнала.

% расчет коэффициентов фильтров:

b = 0.4*sinc(0.4*(-25:25));

brect = b.*rectwin(51)';

bblack = b.*blackman(51)';

% формирование входа:

t=0:9;

x=[exp(t) exp(t) exp(t) exp(t) exp(t) exp(0)];

% фильтрация:

y = filter(b,1,x);

yblack=filter(bblack,1,x);

yrect=filter(brect,1,x);

% вывод в виде графика:

xcoord=0:50;

plot(xcoord,x,xcoord,y,xcoord,yrect,xcoord,yblack)

Как видно в коде, было использовано несколько одинаковых последовательностей экспоненциальной функции. Это делается ввиду быстрого роста функции, во избежание образования большого количества нулей во входной последовательности. x – вектор длиной 51, составленный из пяти экспонент по 10 отсчётов каждая и одного значения экспоненты в нуле.

Переходная характеристика в первой половине отсчётов от «хорошего» фильтра снова оказывается лучше (разница в 2-3%) и, более того, на выходе фильтра, спроектированного с помощью прямоугольного окна, перед каждым ростом экспоненты можно заметить довольно большой спад, которого нет у «хорошего» фильтра. Его величина составляет 4%, если за 100% брать максимальную величину входного сигнала.

Фильтрация функции Хевисайда

Функция Хевисайда может быть определена как:

(Формула 3.2)

Фильтрация функции Хевисайда (с использованием прямоугольного окна и окна Блэкмана)

Рис. 3.12. Фильтрация функции Хевисайда.

b = 0.4*sinc(0.4*(-25:25));

brect = b.*rectwin(51)';

bblack = b.*blackman(51)';

x=ones(51,1); % описание единиц функции при значениях времени больше

% либо равных нулю

% функция ones(m,n) формирует массив единиц размерности m на n

% m – количество строк, n – количество столбцов

y = filter(b,1,x);

yblack=filter(bblack,1,x);

yrect=filter(brect,1,x);

xcoord=0:50;

plot(xcoord,x,xcoord,y,xcoord,yrect,xcoord,yblack)

grid on

На Рисунке 3.12 показана входная функция Хевисайда (синий цвет) и выходы фильтров (красным цветом показан выход фильтра, спроектированного с помощью прямоугольного окна, бирюзовым – выход фильтра, спроектированного с помощью окна Блэкмана).

Если за 100% брать величину входного сигнала, то максимальная разница в переходной характеристике (имеются ввиду первые 25 отсчётов) снова составит 2-3%. Сигнал на выходе «хорошего» фильтра на 2 отсчёта быстрее входит в 10% интервал, чем сигнал на выходе фильтра, спроектированного с помощью прямоугольного окна (при 51 отсчёте эта разница составляет 4%). Если брать 5% интервал, то разница будет уже в 5 отсчётов, что составляет 10% от всех используемых отсчётов.

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

Фильтрация с использованием других окон

Фильтрация функции Хевисайда (с использованием прямоугольного окна и окон Тьюки)

Рис. 3.13. Фильтрация функции Хевисайда. Использовались фильтры, созданные с помощью прямоугольного окна и окон Тьюки.

b = 0.4*sinc(0.4*(-25:25));

brect = b.*rectwin(51)'; % зелёный

bt50 = b.*tukeywin(51,0.5)';

bt75 = b.*tukeywin(51,0.75)';

bt25 = b.*tukeywin(51,0.25)';

x=ones(51,1);

y = filter(b,1,x);

yt50=filter(bt50,1,x);

yt75=filter(bt75,1,x);

yt25=filter(bt25,1,x);

yrect=filter(brect,1,x);

xcoord=0:50;

plot(xcoord,x,xcoord,yrect,xcoord,yt50,xcoord,yt75,xcoord,yt25)

grid

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

Фильтрация функции Хевисайда: переходный процесс (первая половина отсчётов) (с использованием окон Бартлетта-Хеннинга, Бартлетта, Хемминга, Ханна)

Рис. 3.14. Переходный процесс. Использовались Bartlett-Hanning, Bartlett, Hamming, Hann.

b = 0.4*sinc(0.4*(-25:25));

brect = b.*rectwin(51)'; % бирюзовый

bh = b.*barthannwin(51)'; bar = b.*bartlett(51)';

ham = b.*hamming(51)'; han = b.*hann(51)';

x=ones(51,1);

y = filter(b,1,x); yrect=filter(brect,1,x);

ybh=filter(bh,1,x); ybar=filter(bar,1,x);

yham=filter(ham,1,x); yhan=filter(han,1,x);

xcoord=0:50;

xn=0.975*ones(51,1); % для 5% интервала

xv=1.025*ones(51,1); % для 5% интервала

plot(xcoord,x,xcoord,xn,xcoord,xv,xcoord,yrect,xcoord,ybh,xcoord,ybar,

xcoord,yham,xcoord,yhan)

На Рисунке 3.14. изображёна часть переходного процесса. Использовались следующие окна: Bartlett-Hanning (data5 на Рисунке 3.14), Bartlett (data6), Hamming (data7), Hann (data8).

Если единицу брать за 100%, то максимальная разница в первой половине переходной области между всеми сигналами - около 2%.

Фильтрация функции Хевисайда: вторая половина отсчётов (с использованием окон Бартлетта-Хеннинга, Бартлетта, Хемминга, Ханна)

Рис. 3.15. Последние 25 временных отсчётов при фильтровании единичного скачка. Окна Bartlett-Hanning, Bartlett, Hamming, Hann.

Фильтры, полученные с использованием окон Bartlett-Hanning (data5), Hamming (data7), Hann (data8), ведут себя подобно фильтру, полученному применением окна Blackman – они позволяют сигналу раньше войти в 5% интервал (33 отсчёт). Но всё же необходимо заметить, что на 51 отсчёте сигналы data5 и data7 немного отличается от единицы. Сигнал, отфильтрованный с помощью Bartlett, ведёт себя подобно сигналу, отфильтрованному с помощью rectwin – оба входят в 5% интервал на 36 отсчёте. Разница составляет 6% от всех используемых отсчётов. Примечательно, что, из четырёх вышеупомянутых окон, у Bartlett самый высокий уровень первого бокового лепестка, но также нельзя упускать из виду временную область окна Bartlett, существенно отличающуюся от временных областей других окон.

В первой части повествования отмечалось, что треугольное окно и окно Бартлетта имеют похожие графики во временной области, а также одинаковые полосы пропускания, но разные полосы подавления в частотной области. Если взглянуть на график выхода фильтров, спроектированных с помощью этих окон (Рисунок 3.16), то можно увидеть незначительные различия (фильтр, спроектированный с помощью оконной функции bartlett, немного лучше – сигнал на выходе в некоторых местах ближе к единице).

Фильтрация функции Хевисайда: вторая половина отсчётов (с использованием окон Бартлетта и треугольного)

Рис. 3.16. Последние 20 временных отсчётов сигнала на выходе фильтров, спроектированных с помощью bartlett (бирюзовый) и triang (красный).

Оценим фильтрацию с помощью окон Nuttall и Blackman-Harris:

Фильтрация функции Хевисайда: вторая половина отсчётов (с использованием окон Блэкмана-Харриса и Натталла)

Рис. 3.17. Фильтрация с помощью Nuttall (красный) и Blackman-Harris (бирюзовый), вторая половина отсчётов.

b = 0.4*sinc(0.4*(-25:25));

brect = b.*rectwin(51)'; % зелёный

bnut = b.*nuttallwin(51)'; bbh = b.*blackmanharris(51)';

x=ones(51,1);

yrect=filter(brect,1,x);

ynut=filter(bnut,1,x); ybh=filter(bbh,1,x);

xcoord=0:50;

plot(xcoord,x,xcoord,yrect,xcoord,ynut,xcoord,ybh)

grid on

Несмотря на довольно существенную разницу в падении боковых лепестков (Рисунок 1.20), разницы в отфильтрованных сигналах не видно.

Выводы для непараметрических окон (временная область)

Была обнаружена разница между фильтрами, спроектированными с помощью разных типов окон:

Один из самых значительных обнаруженных показателей – поведение переходной полосы (это первые отсчёты, по длине равные половине импульсной характеристики КИХ-фильтра). Максимальное различие между отфильтрованными сигналами – 2-3%.
Вхождение в 5% коридор при установлении постоянного уровня сигнала также может являться показателем окна. Разница между отсчётами при вхождениях сигналов на выходах фильтров в коридор, отнесённая ко всем отсчётам составляет до 10%.
Фильтры, спроектированные с помощью разных окон, по-разному передают входной сигнал, влияя на оценку амплитуды (между отфильтрованными сигналами была обнаружена разница в 1%) и на само представление входной функции после фильтрации (при использовании фильтра, созданного с помощью rectwin, были обнаружены спады перед экспоненциальным ростом функции, причём разница влияния плохого и хорошего окна достигает 4%).

Необходимо понимать, что прямого влияния боковые лепестки на временную область отфильтрованного сигнала не оказывают, так как новый сигнал есть свёртка старого с импульсной характеристикой фильтра. Но, тем не менее, характеристики окна можно определить по максимальному уровню первого бокового лепестка, возможно сделать вывод о том, как будет вести себя фильтр, взвешенный данным окном: как фильтр, импульсная характеристика которого спроектирована с помощью rectwin, или как фильтр, созданный лучшим, по сравнению с прямоугольным, окном.

Рассмотрим уровень первого бокового лепестка окон, принимавших участие в проектировании, заранее разделив окна на подобные по поведению во временной области отфильтрованного сигнала с помощью прямоугольного окна и подобные по поведению отфильтрованного сигнала с помощью окна Блэкмана. Не будем оценивать первый боковой лепесток относительно главного, как это принято в литературе, а оценим его абсолютно, по шкале ординат, в Window Design & Analysis Tool. Это можно с лёгкостью осуществить ввиду нормирования максимального значения временной области (единица).

Окна, подобные прямоугольному по своему поведению, и абсолютный уровень их первого бокового лепестка (данные для 64 отсчётов):

Bartlett (+3,44 дБ)
Triangular, треугольное (+3,55)
Tukey, для параметра 0,5 (+18,36)
Само окно Rectangular, прямоугольное (+22,87)

Окна, подобные окну Блэкмана по своему поведению, и абсолютный уровень их первого бокового лепестка (данные для 64 отсчётов):

Nuttall (-70,65)
Blackman-Harris (-74,87)
Hamming (-16,3)
Hann (-1,53)
Bartlett-Hanning (-5,89)
Само окно Blackman (-29,69)

Итак, у окон, подобных по поведению во временной области хорошему окну Блэкмана, наблюдается отрицательный абсолютный уровень боковых лепестков. У окон, подобных по поведению во временной области прямоугольному окну, наблюдается положительный абсолютный уровень боковых лепестков. Стоит отметить, что уровень главного лепестка (mainlobe, своего рода полоса пропускания окна) варьируется от 27,2 (окно Nuttall) до 36,12 (прямоугольное окно), параметрические окна оценивались с параметром по умолчанию.

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

Окна, влияние которых во временной области не описывается знаком абсолютного уровня первого бокового лепестка

Некоторые окна с параметрами не полностью описываются знаком абсолютного уровня первого бокового лепестка.

Рассмотрим окно Гаусса с разным параметром α (Рисунок 3.18). Чем больше α, тем ниже уровень боковых лепестков АЧХ.

Фильтрация функции Хевисайда: вторая половина отсчётов (с использованием окон Гаусса)

Рис. 3.18. Фильтрация функции Хевисайда с помощью окон Гаусса.

b = 0.4*sinc(0.4*(-25:25));

brect = b.*rectwin(51)'; % зелёный цвет на графике

% далее будет приведён абсолютный уровень боковых лепестков для

% каждого использованного окна (51 отсчёт) и цвет на графике

bg05 = b.*gausswin(51,0.5)'; % +19,76, красный

bg15 = b.*gausswin(51,1.5)'; % +8,84, бирюзовый

bg18 = b.*gausswin(51,1.8)'; % +0.77, фиолетовый

bg19 = b.*gausswin(51,1.9)'; % -3.5, жёлтый

bg25 = b.*gausswin(51,2.5)'; % -16.39, серый

x=ones(51,1);

yrect=filter(brect,1,x);

yg05=filter(bg05,1,x);

yg15=filter(bg15,1,x);

yg18=filter(bg18,1,x);

yg19=filter(bg19,1,x);

yg25=filter(bg25,1,x);

xcoord=0:50;

xn=0.975*ones(51,1);

xv=1.025*ones(51,1);

plot(xcoord,x,xcoord,yrect,xcoord,yg05,xcoord,yg15,xcoord,yg18,xcoord,yg19,xcoord,yg25,xcoord,xn,xcoord,xv); grid on

Из графика (Рисунок 3.18) видно, что чем ниже уровень боковых лепестков применяемого окна, тем больше сигнал на выходе фильтра напоминает сигнал на выходе фильтра, спроектированного с помощью окна Блэкмана. Если абсолютный уровень первого из боковых лепестков не больше 16.01 децибел (α=1), то в случае окна Гаусса отфильтрованный сигнал, если параметром брать вхождение в 5% коридор, напоминает сигнал на выходе фильтра, спроектированного с помощью окна Блэкмана. Окно Гаусса – это, в некотором смысле, исключение из правила. Полезным будет заметить, что максимальная разность в первых 25 отсчётах отфильтрованного сигнала в 2% достигается при параметре α=2.5 и соответствующем абсолютном уровне первого бокового лепестка в -16.39 дБ.

Рассмотрим окно Чебышева. Параметр SidelobeAtten указывает уровень спада боковых лепестков по отношению к главному лепестку окна.

Фильтрация функции Хевисайда: вторая половина отсчётов (с использованием окон Чебышева)

Рис. 3.19. Фильтрация функции Хевисайда с помощью окон Чебышева.
b = 0.4*sinc(0.4*(-25:25));

brect = b.*rectwin(51)'; % зелёный цвет на графике

% далее будет приведён абсолютный уровень боковых лепестков для

% каждого использованного окна (51 отсчёт) и цвет на графике

bc20 = b.*chebwin(51,20)'; % + 5.23, красный

bc25 = b.*chebwin(51,25)'; % +4.91, бирюзовый

bc26 = b.*chebwin(51,26)'; % +4.79 дБ, сиреневый

bc26d5 = b.*chebwin(51,26.5)'; % +4.75 дБ, жёлтый

bc100 = b.*chebwin(51,100)'; % -74.8 дБ, серый

x=ones(51,1);

yrect=filter(brect,1,x);

yc20=filter(bc20,1,x);

yc25=filter(bc25,1,x);

yc26=filter(bc26,1,x);

yc26d5=filter(bc26d5,1,x);

yc100=filter(bc100,1,x);

xcoord=0:50;

xn=0.975*ones(51,1); % для построения 5% интервала

xv=1.025*ones(51,1); % для построения 5% интервала

plot(xcoord,x,xcoord,yrect,xcoord,yc20,xcoord,yc25,xcoord,yc26,xcoord,yc26d5,xcoord,yc100,xcoord,xn,xcoord,xv)

grid

Из графика (Рисунок 3.19) видно, что сигналы, прошедшие через фильтры, взвешенные окнами Чебышева с параметрами SidelobeAtten, меньшими 26.5 даже не входят в 5% интервал. Следовательно, уровень первого из боковых лепестков для 51 отсчёта должен быть меньше либо равен +4.75, что также является исключением: при уровне первого из боковых лепестков +4.75 дБ, созданный фильтр обрабатывает сигнал так, что тот входит в 5% интервал на 3 отсчёта раньше, чем сигнал от фильтра, взвешенного прямоугольным окном! Максимальная разность отфильтрованных сигналов в первых отсчётах будет плохим показателем качества окна – наибольшая разность с сигналом, отфильтрованным b.*rectwin(51)', у сигнала, полученного фильтром, взвешенным окном с параметром SidelobeAtten=20 (меньше, чем у хороших окон Чебышева).

Окно Кайзера также может иметь хорошие характеристики при малом, но положительном, уровне первого бокового лепестка (например, при β=3). Получается, что все параметрические окна, кроме Тьюки, не должны анализироваться по знаку максимального уровня первого бокового лепестка.

Выводы для параметрических окон (временная область)

Параметрические окна, кроме окон Тьюки, не подчиняются тем закономерностям, что были получены для других, непараметрических, окон пакета MATLAB.

Уровни боковых лепестков, при которых окна показывают хороший результат (хорошее вхождение в 5% интервал) при воздействии на функцию Хевисайда:

Окно Гаусса: ≤+16.01 дБ для 51 отсчёта.
Окно Чебышева: ≤+4.75 для 51 отсчёта.
Окно Кайзера: ≤+6.43 для 51 отсчёта.

4. Спектральный анализ

Влияние окна в частотной и временной области: сравнение

Во второй части повествования было сказано, что окна влияют на просачивание спектральных составляющих и на возможности обнаружения двух близких тонов. Просачивание и возможности обнаружения взаимосвязаны: если в другие частотные области будет просачиваться сильный сигнал, то он с лёгкостью может «похоронить», как принято выражаться в английской литературе, слабый сигнал (следовательно, его будет невозможно обнаружить). Просачивание – более важный показатель качества окна, чем вхождение в 5% коридор или поведение в первых отсчётах отфильтрованной функции Хевисайда.

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

Расчёт спектра мощности в MATLAB

Рассчитать спектр мощности в MATLAB гораздо сложнее, чем временную область сигнала на выходе. Спектр мощности в MATLAB чаще всего рассчитывается с участием какого-либо из окон (методы периодограмм и Велча по умолчанию используют прямоугольное окно либо другое, если оно указано в качестве параметра). Но если нужно рассчитать спектр мощности исходного сигнала без умножения на окно, можно и нужно использовать БПФ – один из первых и простейших методов оценки спектра. Будем придерживаться одного и того же метода оценки спектра мощности как для исходного сигнала, так и для отфильтрованного.

Итак, метод оценки спектра мощности сигнала в среде MATLAB:

Fs = 1024; % Частота дискретизации (семплирования) – частота, с

% которой будут браться отсчёты сигнала.

t = 0:1/Fs:1; % Вектор времени, соответствующий частоте Fs.

x = sin(2*pi*t*200)+0.09*sin(2*pi*t*210.5); % Два наложенных друг на друга

% синусоидальных сигнала с частотами 200 Гц и 210.5 Гц и

% соответствующими амплитудами в спектре мощности 0.5 и 0.0045

% (разница более чем в 100 раз).

% Чтобы оптимизировать БПФ, накачивают сигнал достаточным количеством

% нулей так, чтобы его длина была степенью двойки:

nfft = 2^(nextpow2(length(x))); % Нахождение следующей степени двойки

fftx = fft(x,nfft); % Алгоритм БПФ накачивает x некоторым количеством нулей,

% чтобы его длина была степенью двойки, если длина x меньше nfft;

% в случае, если длина x больше nfft, сигнал x обрезается.

% Значения амплитуды БПФ будут симметричны,

% поэтому первые (nfft+1)/2 точек будут уникальными,

% но это число может оказаться дробным.

NumUniquePts = ceil((nfft+1)/2); % ceil округляет аргумент до целого.

fftx = fftx(1:NumUniquePts); % Вектор значений БПФ обрезается,

% остаются только уникальные значения.

mx = abs(fftx); % В вектор mx записываются абсолютные значения fftx.

% Нужно понимать, что в элементах fftx может быть комплексная часть.

mx = mx/length(x); % Масштабирование

mx = mx.^2; % Для оценки мощности амплитуда берётся в квадрате.

% Следующая часть кода, конструкция if, умножает на 2 компоненты

% спектра, не являющиеся компонентами Найквиста:

if rem(nfft, 2)

mx(2:end) = mx(2:end)*2;

else

mx(2:end -1) = mx(2:end -1)*2;

end

f = (0:NumUniquePts-1)*Fs/nfft; % Расчет вектора частот.

plot(f,mx); % Вывод графика спектра мощности.

title('Power Spectrum of a Sine Wave'); % Вывод названия графика.

xlabel('Frequency (Hz)'); % Вывод величины и размерности по x.

ylabel('Power'); % Вывод величины по y.

Спектр мощности, полученный в MATLAB с помощью БПФ

Рис. 3.20. Спектр мощности сигнала sin(2*pi*t*200)+0.09*sin(2*pi*t*210.5) (см. код). Максимальное значение мощности (амплитуда в квадрате) = 0.5.

Соотношение слабого тона и сильного = 0.898%.

Спектр мощности (масштаб)

Рис. 3.21. Соседние тоны: слабый в присутствии сильного, выбран другой масштаб. График тот же, что и на Рисунке 3.20.

Расчёт спектров мощности отфильтрованных с помощью разных окон сигналов

Для расчёта спектра мощности отфильтрованного сигнала нужно пользоваться таким же методом (во-первых, разные методы дают разные результаты, во-вторых, нет необходимости усложнять задачу), как и для расчёта спектра мощности исходного сигнала. Другими словами, сигнал на выходе просто нужно подставить вместо сигнала на входе. Отфильтруем сигнал (Рисунки 3.20 и 3.21) с помощью разных окон, приведём код для получения спектра мощности сигнала на выходе фильтра, взвешенного прямоугольным окном:

Fs = 1024;

t = 0:1/Fs:1;

x = sin(2*pi*t*200)+0.09*sin(2*pi*t*210.5);

b = 0.4*sinc(0.4*(-25:25)); % Расчёт коэффициентов идеального фильтра

brect = b.*rectwin(51)'; % Взвешивание с помощью прямоугольного окна

y=filter(brect,1,x); % Фильтрование

nfft = 2^(nextpow2(length(y)));

ffty = fft(y,nfft);

NumUniquePts = ceil((nfft+1)/2);

ffty = ffty(1:NumUniquePts);

my = abs(ffty);

my = my/length(y);

my = my.^2;

if rem(nfft, 2)

my(2:end) = my(2:end)*2;

else

my(2:end -1) = my(2:end -1)*2;

end

f = (0:NumUniquePts-1)*Fs/nfft;

plot(f,my);

title('Power Spectrum (rectwin)');

xlabel('Frequency (Hz)');

ylabel('Power');

grid

Любое окно вносит коррективы в амплитуду отсчётов спектра. Далее ко всем спектральным характеристикам прилагается информация о максимальной мощности в спектре на выходе и о процентном соотношении значений мощностей слабого тона и сильного (чем это отношение больше, тем лучше). Интервал в 10.5 дБ между тонами для нашего примера с синусоидами довольно большой, поэтому структура боковых лепестков большого сигнала даже для прямоугольного окна не может поглотить главный лепесток слабого сигнала, из-за чего можно численно сравнить влияние окон на спектр мощности сигнала.