» % диаграмма нулей и полюсов
Для вычисления выходного сигналаy линейной дискретной системы с переджаточной функцией 
при произвольном входном сигнале x может использоваться команда
y = filter(num, den, x, zi),
при этом ненулевые начальные условия задаются в векторе zi. Количество элементов в этом векторе определяется как max(length(num), length(den)-1). В частности, для системы второго порядка с N = M =2 этот вектор равен
zi = [b1*x[-1]+b2*x[-2]-a1*y[-1]-a2*y[-2], b1*x[-1]-a2*y[-1]].
» y=filter(num, den, x, zi);
» stem(n, y), xlabel(‘n’), ylabel(‘y(n)’)
Функция tf2zp() выполняет преобразование коэффициентов полиномов передаточной функции дискретной системы в нули и полюсы передаточной функции, функция zp2tf() осуществляет обратное преобразование.
Частотные характеристики дискретных систем
Расчет частотной характеристики (H(e jω ) дискретной системы выполняется с помощью одной из функций freqz(…). Функция freqz(a, b) вычисляет частотную характеристику и выводит в текущее графическое окно амплитудную частотную (АЧХ) и фазовую частотную (ФЧХ) характеристики.
Функция [h, f]=freqz(num, den, n, Fs) вычисляет n комплексных отсчетов частотной характеристики по векторам num, den при частоте дискретизации Fs. Возвращаемый параметр f содержит отсчеты частоты в герцах, равномерно распределенные в интервале от 0 до Fs/2, количество элементов h и f равно n.
» ylabel(‘Модуль (дБ)’), xlabel(‘Частота, Гц’), grid
Функция Q = unwrap(P) корректирует фазо-частотную характеристику системы при переходе через значения π, устраняя разрывы функции.
Пример. Производится расчет фильтра нижних частот Чебышева 1 типа и строятся его частотные характеристики
Fs=1000; % Частота отсчетов
% Граничные частоты полосы пропускания и задерживания фильтра
% Максимально допустимые пульсации
% в полосе пропускания и задерживания
% Расчет порядка и передаточной функции ФНЧ
[n,Wn]=cheb1ord(Wp,Ws, Rp, Rs);
% Вычисление частотной характеристики
figure(1), subplot(3,1,1), plot(f, abs(h)), title(‘ АЧХ фильтра’)
subplot(3,1,2), plot(f, angle(h)), title(‘ФЧХ фильтра’)
% Построение скоректированой ФЧХ
subplot(3,1,3), plot(f, unwrap(angle(h))), xlabel(‘Частота, Гц’)
Передаточная функция и частотная характеристика линейной дискретной системы часто выражается через нули и полюса. Их принято представлять на комплексной плоскости. Для этого используют функции zplane(z,p) и zplane(b,a).
Вычисление оконных функций
Оконные функции используются для ограничения задания сигнала во временной области или для ограничения частотного спектра сигнала. Важными свойствами окон являются их амплитудные спектры. Фазовые характеристики окон имеют линейный характер и поэтому, как правило, не рассматриваются.
Примеры некоторых оконных функций. Функция w=hanning(n) возвращает n – точечное окно Хэннинга в векторе-столбце w. Коэффициенты окна Хэннинга вычисляются по выражению
.
Аналогично функция 
вычисляет n – точечное окно Блэкмана.
Среди 14 оконных функция Matlab имеются окна с регулируемыми параметрами. Примером является окно Кайзера. Синтаксис соответствующей функции Matlab
Здесь параметр beta определяет минимальное затухание боковых лепестков спектра окна.
Функция w=window(fhandle, n) возвращает n – точечное окно типа, задаваемого параметром fhandle.
Пример вывода трех окон и их амплитудных спектров с помощью специальной функции визуализации окон wvtool(w1,w2,…, wn)
Вычисления с комплексными числами…………………………………………………. 5
Получение справочных сведений в MATLAB…………………………………………. 6
Поэлементные операции с массивами…………………………………………………..12
Программирование в MATLAB…………………………………………………………13
Построение графиков MATLAB………………………………………………………….17
Анализ линейных непрерывных систем………………………………………………. 27
Частотные характеристики и их графики………………………………………………..28
Анализ линейных дискретных систем…………………………………………………. 32
Частотные характеристики дискретных систем………………………………………….35
Вычисление оконных функций…………………………………………………………. 37
Чтобы распечатать файл, скачайте его (в формате Word).
MATLAB\MATLAB
h = freqs(b,a,w)
[h,w] = freqs(b,a)
[h,w] = freqs(b,a,f)
freqs(b,a)
Функция freqs рассчитывает комплексную частотную характеристику H ( j w ) аналоговой системы, заданной своей функцией передачи в s-области (преобразованием Лапласа):
Исходными данными для расчета служат векторы b и a коэффициентов числителя и знаменателя функции передачи.
Автоматически выбирает для расчета частотной характеристики h набор из 200 частотных точек, возвращая их значения во втором выходном параметре w.
Функция freqs, вызванная без указания выходных параметров, строит в текущем графическом окне графики амплитудно-частотной и фазочастотной характеристик.
Функция freqs работает только с вещественными системами и положительными частотами.
Построим графики частотных характеристик аналоговой системы, имеющей следующую функцию передачи:
.
a = [1 0.4 1];
b = [0.2 0.3 1];
w = logspace(-1,1);
freqs(b,a,w)
Графики можно получить также следующим образом:
h = freqs(b,a,w);
mag = abs(h);
phase = angle(h);
subplot(2,1,1), loglog(w,mag)
subplot(2,1,2), semilogx(w,phase)
Чтобы отобразить частоту в герцах, ФЧХ в градусах, а АЧХ в децибелах, можно поступить так:
f = w/(2*pi);
mag = 20*log10(mag);
phase = phase*180/pi;
Функция freqs рассчитывает значения полиномов в заданных частотных точках, а затем делит значения полинома числителя на значения полинома знаменателя:
s = i*w;
h = polyval(b,s)./polyval(a,s);
MATLAB 8.0 (R2012b): создание, обработка и фильтрация сигналов, Signal Processing Toolbox
Введение
Окно справки по пакету расширения Signal Processing Toolbox системы MATLAB 8.0 представлено на рис. 1 на фоне рабочего окна самой системы с открытой вкладкой каталога пакетов расширения APPS. Одна из кнопок в панели инструментов Signal Analysis дает доступ к браузеру сигналов, фильтров и спектров, показанному в правой части окна справки. В левой части этого окна указаны наименования разделов пакета Signal Processing Toolbox.
Рис. 1. Окно справки по пакету расширения Signal Processing Toolbox
Как видно на рис. 1, пакет Signal Processing Toolbox состоит из следующих разделов:
Создание сигналов
Многие сигналы представлены как функции времени s(t), параметры которой можно изменять с помощью модуляции того или иного вида. Модуляцией называют процесс изменения какого-либо параметра (амплитуды, частоты, фазы и т. д.) по определенному закону, в результате чего сигнал становится переносчиком информации.
MATLAB 8.0 со своими встроенными средствами позволяет создавать множество сигналов. Например, простейшим является синусоидальный сигнал:
где A — амплитуда; f — частота; j — фаза сигнала. Задав эти параметры ( t — как вектор отсчетов сигнала, число элементов которого определяет число отсчетов сигнала):
можно легко построить график сигнала s ( t ):
Он показан на рис. 2 в графическом окне системы MATLAB. Такое построение осуществляется в командной строке, на что указывает приглашение к вводу >>. В командном окне также показано окно с краткими данными о системе MATLAB 8.0 (R2012b). Согласно этим данным система выпущена на рынок в 2012 году.
Рис. 2. Окно командного режима
В дальнейшем, учитывая множество сигналов, мы будем представлять их по четыре в каждом графическом окне, используя для этого функцию его разбиения на четыре под-окна subplot.
Рис. 3. Окно с подокнами меандра, апериодического треугольного импульса, пилообразного периодического импульса и периодического треугольного импульса
Представленная программа (она задается в редакторе) дает примеры создания и графической визуализации четырех типов простых сигналов (меандра, одиночного треугольного импульса, пилообразного импульса и симметричного треугольного импульса) (рис. 3):
Рис. 4. Временные диаграммы четырех сигналов, относящихся к функциям Гаусса
Важное значение имеет функция sinc (или sin(π t )/π t при t ≠ 0 и 1 при t = 0). Функция sinc( t ) представляет обратное преобразование Фурье для прямоугольного импульса с высотой 1 и шириной 2π:
Кроме того, эту функцию можно использовать как базисную для восстановления любого сигнала g ( t ) по его отсчетам, если спектр сигнала ограничен условием –p :
формирует выборку (дискретные значения) косинусоидального сигнала с частотой от f 0 в начальный момент времени t до f 1 в конечный момент времени t 1. Звук такого сигнала напоминает визг, откуда и его название (chirp). По умолчанию t = 0, f 0 = 0 и f 1 = 100. Необязательный параметр phi (по умолчанию 0) задает начальную фазу сигнала. Другой необязательный параметр — method — задает закон изменения частоты: linear — линейный закон (по умолчанию), quadratic — квадратичный и logarithmic — логарифмический.
Функция strips позволяет детально (по частям с длительностью 0,25 с) рассмотреть первые два сложных сигнала, например частотно-модулированного от функции vco. В разделе модуляция/демодуляция есть пример с GUI (рис. 7), иллюстрирующий такие сигналы с различными видами модуляции.
Рис. 7. Визуализация модулированных сигналов
Свертка и корреляция сигналов
При записи этого выражения учтено, что нумерация индексов массивов в MATLAB идет с единицы (1). Свертка реализуется функцией conv(a,b):
Операцию свертки часто используют для вычисления сигнала на выходе линейной системы y по сигналу на входе x при известной импульсной характеристике системы h :
Эту операцию можно использовать для осуществления простейшей фильтрации сигнала:
Для двух векторов x и y с длиной m и n определена операция свертки:
Рис. 8. Линейная и циклическая свертка
Обратная свертке функция — это [q,r] = deconv(z,x). Она фактически определяет импульсную характеристику фильтра (рис. 8):
Рис. 9. Очистка от шума линейно-нарастающего сигнала
Свертку часто применяют для очистки сигналов от шума (рис. 9):
Для двумерных массивов также существует функция свертки: Z = conv2(X,Y) и Z = conv2(X,Y,‘option’). Возможна и многомерная свертка — функция convn. Новая функция mscohere строит график зависимости квадрата модуля функции когерентности от частоты (рис. 10):
Рис. 10. Зависимость квадрата модуля функции когерентности от частоты
Следующая программа строит отсчеты двух сигналов с задержкой на три отсчета — треугольного сигнала и шума (рис. 11):
Рис. 11. Два коррелированных сигнала — треугольный и шума
Кросс-корреляцию этих сигналов обеспечивает программа (рис. 12):
Рис. 12. Кросс-корреляция сигналов
Преобразование сигналов
В MATLAB всегда много внимания уделялось различным преобразованиям сигнала [2–5], в частности спектральным. Спектр дискретного сигнала является периодическим, и прямое дискретное преобразование Фурье (ДПФ или Discrete Fourier Transform, DFT) определяется выражением:
Для предотвращения растекания (размазывания) спектра дискретных сигналов часто используются окна. Для этого достаточно в формуле прямого ДПФ под знаком суммы ввести еще один множитель — W ( k ). Соответственно, обратное дискретное преобразование Фурье задается выражением:
Частотным спектром случайного процесса является преобразование Фурье от корреляционной функции случайного процесса Rx :
В радиоэлектронике особый интерес представляет спектральная оценка сильно зашумленных сигналов. Для таких сигналов применяются два подхода: непараметрический — использующий только информацию, извлеченную из сигнала (реализован в методах периодограмм и Уэлча), и параметрический — предполагающий наличие некоторой статистической модели сигналов, параметры которой подлежат определению. Реализовано восемь классов алгоритмовспектрального анализа : Periodogram, Welch, MTM (Thomson multitaper method), Burg, Covariance, Modified Covariance, Yule-Walker, MUSIC (Multiple Signal Classification) и Eigenvector. Их подробное описание дано в [4–6].
Периодограммы, спектрограммы и их применение
Обычный спектр строится методом быстрого преобразования Фурье (БПФ, FFT) часто с применением временного окна, предотвращающего разрывы сигнала на концах интервала анализа спектра. Например, так реализованы периодограммы. Для построения спектрограмм используется разбивка интервала анализа на короткие окна. Короткое окно пробегает общий интервал анализа, и в каждом частичном интервале строится свой спектр. Их наложение дает спектрограмму, определенную в пространстве «уровень – частота – время» (на плоскости уровень представляется цветом), тогда как обычный спектр определен в плоскости «уровень – частота».
Рис. 13. Построение спектрограмм сигналов, модулированных по различным законам: логарифмическому, линейному, треугольному и синусоидальному
На рис. 13 показано построение четырех спектрограмм для сигналов с различными законами частотной модуляции. Нетрудно увидеть, что во всех случаях закон модуляции отчетливо распознается и позволяет судить об области частот, в которой действует сигнал (chirp или vco). Этот рисунок строит следующая программа:
Рис. 14. Сравнение периодограммы и спектрограммы сигнала sinc(t)
Из сказанного может сложиться неверное представление о явных преимуществах спектрограмм по сравнению с периодограммами (функция psd). То, что это далеко не так, показывает программа для сигнала sinc (рис. 14):
Эта программа строит периодограмму и спектрограмму функции sinc(t). Теперь беспомощной оказывается спектрограмма, по которой ничего нельзя сказать о спектре сигнала и области его частот. А периодограмма более информативна: она указывает на вид спектра и область занимаемых им частот. В частности, хорошо видно постоянство спектра в начальной области частот, присущее этой функции. Правда, вид спектрограммы сильно зависит от типа короткого окна: в данном случае задано окно Блэкмана (Blackman). Вид периодограммы также зависит от выбора окна, но глобального.
Чем сложнее сигнал, тем более детальной и эффектной оказывается спектрограмма. На рис. 15 показан пример модуляции/демодуляции с построением спектрограммы сложного звукового сигнала. Кстати, оригинальный сигнал и сигнал, прошедший модуляцию/демодуляцию, можно воспроизвести на компьютере, оборудованном звуковой картой и акустической системой. Это можно сделать при различных видах модуляции.
Рис. 15. Спектрограмма сложного акустического сигнала с амплитудной модуляцией
Похожие на спектрограммы картинки дают вейвлетограммы и скайлеграммы, получаемые при вейвлет-анализе сигналов [7]. Порою они более информативны. Но вейвлеты в Signal Ptocessing Toolbox не реализованы: они описаны и используются в отдельном пакете расширения Wavelet Toolbox.
Для выполнения дискретного преобразования Фурье различными методами служит GUI-окно Discrete Fourier Transform (рис. 16).
Рис. 16. Спектр пилообразного сигнала в окне дискретного Фурье-преобразования
Оно позволяет задать один из трех видов сигнала (синусоида, меандр и пилообразный) и его периодограмму (спектр) при одном из шести видов окон. Как сигнал, так и окно можно загружать извне из файла или рабочего пространства MATLAB. Нетрудно убедиться в большом влиянии на вид спектра выбранного окна.
Оконные функции и браузер окон
Учитывая важную роль окон, в Signal Processing Toolbox входит 21 N ‑точечное окно ( N — целое число), называемое по фамилии предложивших окно ученых, например Hamming, Blackman, Bartlett, Chebyshev, Taylor, Kaiser и т. д. Для просмотра временных и амплитудно-частотных характеристик всех окон есть соответствующие функции, но удобно пользоваться GUI-браузером окон wvtool.
позволяет строить четыре типа 64‑точечных окон — прямоугольное, Хемминга, Ханна и Гаусса. Можно задать окна и отдельными командами.
Рис. 17. Сравнение в браузере трех окон Кайзера с параметром бетта 1,5, Блэкмана и Блэкмана-Харисса
Это, а также сравнительное построение трех других типов окон обеспечивает следующая программа (рис. 17):
Рис. 18. Конструктор окон
Существует также конструктор-анализатор окон с GUI-интерфейсом, окно которого (рис. 18) открывается командой:
По умолчанию в нем открывается окно объекта sigwin, но можно открыть и другие окна (в том числе окно пользователя) из списка Current Window Information. Открытые окна появляются в третьем нижнем окне. Окна можно скопировать, добавить в список, установить в рабочее пространство MATLAB или стереть.
Изменение числа отсчетов и интерполяция сигналов
Изменение числа отсчетов широко используется в технике цифровой обработки сигналов. Для уменьшения числа отсчетов применяются операция децимации и функция decimate (рис. 19):
Рис. 19. Уменьшение частоты дискретизации сигнала (децимация)
Для увеличения числа отсчетов исходный сигнал интерполируют (функция interp), а затем нужное число отсчетов сигнала берут из кривой интерполированного сигнала (рис. 20):
Рис. 20. Интерполяция сигнала и увеличение частоты его дискретизации
Аналоговые и цифровые фильтры
Фильтры имеют особое значение при обработке сигналов. С их помощью осуществляется очистка сигналов от шума или реализуются избирательные свойства систем. Будем считать, что читатель знаком с теорией фильтров.
Фильтрующие цепи обычно задаются своей операторной передаточной характеристикой:
Имея векторы коэффициентов полиномов a ( s ) и b ( s ), с помощью функции freqs можно построить АЧХ и ФЧХ фильтрующей цепи в логарифмическом масштабе (рис. 21):
Рис. 21. Построение логарифмической АЧХ и ФЧХ системы по ее операторной передаточной характеристике
В Signal Processing Toolboox входит множество функций по расчету и проектированию различных фильтров — нижних, верхних частот и полосовых. Порою для получения важных характеристик фильтров достаточно задать нужную строку программного кода. Например, построение характеристик аналогового фильтра Бесселя (рис. 22) реализуется следующим программным фрагментом:
Для просмотра практически всех характеристик фильтров можно использовать визуализатор фильтров fvtool. Приведенная ниже программа дает характеристики фильтра Бесселя нижних частот (рис. 23):
Рис. 23. Характеристики аналогового эллиптического фильтра нижних частот
Открытое меню Analysis дает представление обо всех доступных видах анализа фильтров. Среди них построение АЧХ и ФЧХ фильтра, импульсная и переходная характеристики, групповая задержка и т. д.
Конструирование двух цифровых фильтров НЧ (FIR и Баттерворта) со сравнительным построением их характеристик в окне fvtool (рис. 24) обеспечивает следующая программа:
Рис. 24. АЧХ двух цифровых фильтров для сравнения
Программа строит и другие характеристики фильтров, например переходные и импульсные. Напомним, что переходная характеристика является реакцией системы (фильтра) на единичный скачок, а импульсная — на импульс единичной площади с длительностью, стремящейся к нулю (рис. 25).
Рис. 25. Импульсные характеристики двух цифровых фильтров
Интерактивный конструктор-анализатор фильтров fdatool
В Signal Proccesing Toolbox входит интерактивный конструктор-анализатор фильтров с GUI-интерфейсом. Он позволяет без какого-либо программирования анализировать и проектировать семь основных типов фильтров (таблица).