ыс. А. Б. Сергиенко минобрнауки россии санктПетербургский государственный электротехнический университет лэти им. В. И. Ульянова (Ленина) А. В. Петров а. Б. Сергиенко цифровая обработка сигналов лабораторный практикум
 Скачать 1.73 Mb.
|
|
3.3. Индивидуальное задание В данной лабораторной работе используется сигнал из индивидуального задания для лабораторной работы № 1. Дополнительные индивидуальные па- раметры отсутствуют. 3.4. Указания к выполнению работы 1. Подготовка к началу работы. Запустите MATLAB и сделайте теку- щей папку вашей бригады на сервере дисплейного класса. Создайте в редакторе MATLAB новый файл MATLAB-программы и скопируйте в него часть кода из программы, созданной при выполнении ла- бораторной работы № 1 или № 2. Необходимо выбрать фрагмент кода, с по- мощью которого формируется исходный дискретный сигнал и строится его график. Сохраните новую программу с соответствующим именем и добейтесь ее стабильной работы (программа должна формировать вектор отсчетов дис- кретного сигнала и строить его график). Последующие пункты работы выполняются путем дополнения создан- ной MATLAB-программы. 2. Расчет ДПФ. Вычислите ДПФ для дискретного сигнала из индивиду- ального задания. Используя функцию stem, постройте графики модуля и фа- зы спектральных отсчетов в одном графическом окне друг под другом (с по- мощью функции subplot). 39 Справка Для вычисления ДПФ в MATLAB имеется функция fft (она вычисляет ДПФ с использованием быстрых алгоритмов — FFT, Fast Fourier Transform). В простейшем случае она вызывается без дополнительных параметров: y = fft(x) . Здесь x — вектор отсчетов сигнала; y — вектор результатов вычисления ДПФ. Замечание Чтобы улучшить внешний вид частотных графиков, целесообразно изменить принятые по умолчанию пропорции графического окна, растянув его в ширину. 3. Оценка ширины спектра сигнала. Используя обратное ДПФ, опре- делите минимальное число низкочастотных гармонических составляющих сигнала, содержащих не менее 90 % его энергии. Примерная последователь- ность действий: 1. Рассчитайте энергию исходного сигнала: 2 0 ( ) k E x k = ∑ (3.6) 2. Создайте копию вектора результатов ДПФ под новым именем. 3. Обнулите в созданной копии результатов ДПФ те элементы, которые соответствуют гармоникам с номерами, превышающими некоторое порого- вое значение N max . Исходное значение N max следует принять равным 0. Замечание При определении диапазона номеров обнуляемых элементов ДПФ следует принять во внимание следующее: • нумерация элементов массивов в MATLAB начинается с единицы; • первое значение в массиве результатов ДПФ соответствует постоянной со- ставляющей сигнала ( (0) X ɺ ); • для получения корректного спектра вещественного сигнала в векторе ре- зультатов ДПФ должны сохраняться пары значений, соответствующие парам одинаковых по модулю положительных и отрицательных частот; • вторая половина вектора результатов ДПФ с учетом периодичности спек- тров дискретных сигналов соответствует отрицательным частотам. В каче- стве примера в табл. 3.1 приведены разные способы нумерации элементов 8- точечного ДПФ. Двойной рамкой и полужирным шрифтом в таблице выделен диапазон элементов вектора, который необходимо обнулить в случае N max = 2. 40 Таблица 3.1 Номер элемента вектора MATLAB 1 2 3 4 5 6 7 8 Номер n из теоретической формулы ДПФ 0 1 2 3 4 5 6 7 Нумерация с использованием отрицательных частот 0 1 2 3 ±4 −3 −2 −1 4. Вычислите обратное ДПФ и постройте график получившегося сигна- ла функцией stem. Справка Вычисление обратного ДПФ производится с помощью функции ifft: x = ifft(y); 5. Вычислите энергию получившегося сигнала по (3.6). Если она мень- ше, чем 90 % энергии исходного сигнала E 0 , увеличьте N max на единицу и повторите шаги 3–5. Поиск значения N max можно реализовать вручную (задавая новое значе- ние, запуская программу и анализируя результаты) либо написать MATLAB- программу, которая будет это делать автоматически. В результате выполнения данного пункта работы должно быть получено значение N max , при котором сигнал после обратного БПФ содержит не менее 90 % энергии исходного сигнала. Кроме того, должен быть построен график указанного сигнала вместе с исходным в общих координатных осях. 4. Дополнение сигнала нулями. Скопируйте исходный дискретный сигнал в переменную с новым именем и добавьте к концу этой копии сигнала нулевые отсчеты в количестве, равном длине сигнала (длина вектора должна, таким образом, увеличиться в 2 раза). Вычислите ДПФ для дополненного ну- лями сигнала, постройте (с помощью функции stem) графики модуля и фазы спектральных отсчетов в одном графическом окне друг под другом (с помо- щью функции subplot). Замечание Чтобы улучшить внешний вид частотных графиков, целесообразно изменить принятые по умолчанию пропорции графического окна, растянув его в ширину. 5. Измерение скорости расчетов при вычислении ДПФ непосредст- венно по теоретической формуле. Функция fft гибко выбирает алгоритм вычисления в зависимости от длины преобразуемого сигнала, поэтому изме- рить с ее помощью скорость вычислений по прямой формуле не удастся. Для реализации такого измерения придется вычислять ДПФ как произвольное 41 линейное преобразование — путем перемножения вектора сигнала и квад- ратной матрицы преобразования. Справка Матрица преобразования для ДПФ рассчитывается функцией dftmtx (N — размер вычисляемой матрицы ДПФ): D = dftmtx(N); Для измерения скорости вычислений необходимо реализовать следую- щий код (предполагается, что вектор сигнала x — строка). Для измерения времени выполнения фрагмента кода используется пара функций tic и toc (первая запускает таймер, вторая выводит результат), вычисления повторяют- ся в цикле много раз, чтобы время выполнения кода было достаточно велико. Внимание! Число повторений цикла (в приведенном примере — 1000) необходимо по- добрать экспериментально, чтобы время расчета при N = 1024 составляло примерно 1 секунду. Это время ощутимо варьируется от запуска к запуску, так что стремиться к большой точности здесь не следует. N = 1024; % размер ДПФ % дополнение сигнала нулями до длины N x1 = [x zeros(1, N-length(x))]; D = dftmtx(N); % матрица ДПФ y = zeros(1, N); % массив для результатов ДПФ tic % старт таймера for k = 1:1000 % цикл для измерения времени y = x1 * D; % вычисление ДПФ по прямой формуле end toc % отображение измеренного времени Запишите подобранное число повторений цикла K и занесите измерен- ное время выполнения кода при N = 1024 в табл. 3.2. Таблица 3.2 N Параметр 64 128 256 512 1024 2048 4096 8192 Время, с Время одно- кратного вычис- ления ДПФ, мкс Повторите эксперимент для остальных значений N, указанных в табли- це, не меняя подобранное число повторений цикла K. 42 Внимание! Все измерения должны осуществляться на одном и том же компьютере. Последняя строка таблицы заполняется при подготовке отчета путем де- ления измеренных значений времени на использованное число повторений цикла K. Внимание! Время работы кода растет пропорционально N 2 , поэтому для последнего зна- чения 8192 оно будет составлять около 1 минуты. 6. Измерение скорости расчетов при вычислении ДПФ с использо- ванием быстрого алгоритма. Функция fft гибко выбирает алгоритм вы- числения в зависимости от длины преобразуемого сигнала, и при длине, рав- ной степени двойки, будет использован широко известный алгоритм Кули— Тьюки, число операций в котором пропорционально N log 2 (N). Для измерения скорости вычислений необходимо реализовать следую- щий код. Его общая структура аналогична коду, приведенному выше, однако в явном виде дополнять сигнал нулями не нужно — функции fft можно пе- редать дополнительный входной параметр, указывающий размерность вы- числяемого ДПФ: y = fft(x, N). Внимание! Число повторений цикла (в приведенном примере — 100 000) необходимо подобрать экспериментально, чтобы время расчета при N = 1024 составляло примерно 1 секунду. Это число должно оказаться примерно в 100 раз больше, чем в предыдущем случае. N = 1024; % размер ДПФ y = zeros(1, N); % массив для результатов ДПФ tic % старт таймера for k = 1:100000 % цикл для измерения времени y = fft(x, N); % вычисление БПФ end toc % отображение измеренного времени Запишите подобранное число повторений цикла K и занесите измерен- ное время выполнения кода при N = 1024 в таблицу, по форме повторяющую табл. 3.2. 43 Повторите эксперимент для остальных значений N, указанных в табли- це, не меняя подобранное число повторений цикла K. Внимание! Все измерения должны осуществляться на том же компьютере, на котором производились измерения при вычислении ДПФ непосредственно по теоре- тической формуле. Последняя строка таблицы заполняется при подготовке отчета путем де- ления измеренных значений времени на использованное число повторений цикла K. 7. Подготовка материалов для отчета. В завершение работы необхо- димо скопировать в документ Microsoft Word (в качестве заготовки для отче- та) созданный программный код и все полученные графики: • исходный дискретный сигнал; • модуль и фазу результатов ДПФ; • дискретный сигнал, восстановленный после обнуления части гармо- ник, и исходный сигнал, построенные в общих координатных осях; • модуль и фазу результатов ДПФ при дополнении сигнала нулями. 3.5. Содержание отчета Отчет должен содержать созданный в процессе работы программный код, оформленный в виде законченного документа (с заголовками разделов, формулами и комментариями к коду). Полученные в ходе работы графики и заполненные таблицы размещаются в соответствующих разделах отчета. Используя значение N max , необходимо рассчитать ширину спектра сиг- нала в герцах. Кроме того, в соответствующем разделе должны быть построены графи- ки затрат времени на однократное вычисление ДПФ прямым и быстрым ме- тодами. Для получения этих данных следует разделить занесенные в таблицы значения времени на использованное количество повторений цикла. Вдоль горизонтальной оси следует откладывать log 2 (N). Вертикальная ось также должна быть представлена в логарифмическом масштабе. Совместно с этими экспериментальными графиками следует построить аппроксимирующие зависимости: • для прямого расчета ДПФ: t = k 1 N 2 ; • для быстрого расчета ДПФ: t = k 2 N log 2 (N). 44 Коэффициенты k 1 и k 2 необходимо подобрать самостоятельно так, что- бы аппроксимирующая зависимость при значениях N = 1024 и 8192 была как можно ближе к результатам измерений. Значения k 1 и k 2 должны быть при- ведены в отчете. В конце отчета необходимо привести выводы по результатам работы. 3.6. Контрольные вопросы 1. Запишите матрицу ДПФ для N = 4. 2. Как вычислить энергию сигнала через его ДПФ? 3. Является ли монотонной зависимость энергии сигнала после обнуле- ния части спектра и обратного ДПФ от числа использованных гармоник? От- вет обосновать. 4. Можно ли найти такой сигнал {x(k)} длиной N отсчетов, чтобы его ДПФ совпадало с самим сигналом, т. е. чтобы для всех n выполнялось равен- ство ( ) ( ) X n x n = ɺ ? 5. Длину дискретного сигнала увеличили в 2 раза путем дублирования каждого отсчета ({x(0), x(0), x(1), x(1), …, x(N − 1), x(N − 1)}). Как изменятся результаты ДПФ? 6. Сигнал представляет собой сумму двух синусоид с частотами 300 и 320 Гц. Частота дискретизации равна 10 кГц, длина сигнала 500 отсчетов. Как примерно будет выглядеть модуль ДПФ этого сигнала? 7. Сигнал представляет собой сумму двух синусоид с частотами 400 и 450 Гц. Частота дискретизации равна 5 кГц, длина сигнала 200 отсчетов. Как примерно будет выглядеть модуль ДПФ этого сигнала? 8. Сигнал представляет собой сумму двух синусоид с частотами 500 и 525 Гц. Частота дискретизации равна 10 кГц, к сигналу длиной 400 отсчетов добавлено столько же нулевых значений. Как примерно будет выглядеть мо- дуль ДПФ этого сигнала? 9. Последовательность отсчетов {x(k)} длиной N (k = 0, 1, …, N − 1) подвергли прямому ДПФ. К полученному результату еще раз применили прямое ДПФ. Чему будет равен результат? (Выразить его через {x(k)}.) 10. Последовательность отсчетов {x(k)} длиной N (k = 0, 1, …, N − 1) подвергли обратному ДПФ. К полученному результату еще раз применили обратное ДПФ. Чему будет равен результат? (Выразить его через {x(k)}.) 45 11. Отсчеты последовательности {x(k)} являются чисто мнимыми. Каки- ми свойствами благодаря этому будет обладать ее ДПФ? 12. Придумайте максимально эффективную (по числу арифметических операций) схему реализации ДПФ для N = 3. Сколько вещественных сложе- ний/вычитаний и умножений она требует? (Считать, что вычитание — от- дельная операция, по сложности эквивалентная сложению, так что операции умножения на − 1 не требуются.) 13. Как с помощью ДПФ можно получить N равномерно расположенных отсчетов спектра ( ( ) 2 X j n N π ɺ , n = 0, 1, …, N − 1) последовательности ко- нечной длины, состоящей из M > N отсчетов? 14. Запишите матрицу ДПФ для N = 6. 15. К спектральному отсчету ( 2) X N ɺ прибавили единицу. Как изменит- ся последовательность отсчетов {x(k)}, которой это ДПФ соответствует? 16. Экспериментальные и аппроксимирующие зависимости временных затрат на однократное вычисление ДПФ прямым и быстрым методами рас- ходятся при низких значениях N. Как это можно объяснить? 46 4. ПРОЕКТИРОВАНИЕ ДИСКРЕТНЫХ ФИЛЬТРОВ 4.1. Цели работы Синтез дискретного ФНЧ, удовлетворяющего заданным требованиям, путем использования нескольких алгоритмов: • билинейного преобразования и четырех стандартных аналоговых про- тотипов; • синтеза с использованием оконных функций; • минимаксного синтеза; • минимизации среднеквадратической ошибки. 4.2. Теоретические сведения Под проектированием (или синтезом) цифрового фильтра понимается выбор таких наборов коэффициентов {a i } и {b i }, при которых характеристи- ки получающегося фильтра удовлетворяют заданным требованиям. При синтезе ФНЧ набор задаваемых параметров может выглядеть сле- дующим образом (рис. 4.1): • F д — частота дискретизации; • F pass — граница полосы пропускания; • F stop — граница полосы задерживания; • A pass — допустимая неравномерность АЧХ в полосе пропускания; • A stop — требуемое подавление сигнала в полосе задерживания. 0 f | K( ) f | F pass A pass A stop F stop F д /2 Рис. 4.1. Параметры, требуемые для синтеза ФНЧ Серые области на рисунке демонстрируют допуски, в которые должна укладываться АЧХ фильтра в полосах пропускания и задерживания. В пере- 47 ходной зоне (диапазон частот F pass …F stop ) никаких требований к АЧХ рас- считываемого фильтра не предъявляется. Номинальное значение коэффициента передачи фильтра в полосе про- пускания, как правило, равно единице (0 дБ). В зависимости от способа рас- чета фильтра отклонения коэффициента передачи в полосе пропускания от единицы могут быть как односторонними, так и двусторонними. Методы синтеза цифровых фильтров можно классифицировать по раз- личным признакам: • по типу получаемого фильтра: ‒ методы синтеза рекурсивных фильтров; ‒ методы синтеза нерекурсивных фильтров; • по наличию аналогового прототипа: ‒ методы синтеза с использованием аналогового прототипа (метод би- линейного z-преобразования); ‒ прямые (без использования аналогового прототипа) методы синтеза: оптимальные, в которых численными итерационными метода- ми ищется минимум заданной функции качества (минимаксный син- тез и минимизация квадратической ошибки); субоптимальные, не дающие оптимального решения, но позво- ляющие значительно упростить вычисления по сравнению с опти- мальными методами (синтез с использованием окон). Метод билинейного z-преобразования позволяет синтезировать рекур- сивный дискретный фильтр по частотной характеристике аналогового прото- типа. При синтезе дискретного фильтра по аналоговому прототипу необхо- димо реализовать переход из p-области в z-область, т. е. преобразовать функ- цию передачи аналогового фильтра H(p) в функцию передачи дискретного фильтра H(z). Получающийся дискретный фильтр не может быть полностью идентичен аналоговому по своим характеристикам — хотя бы потому, что частотные характеристики дискретного фильтра являются периодическими. Можно говорить только об определенном соответствии характеристик ана- логового и дискретного фильтров. Функция передачи аналоговой цепи с сосредоточенными параметрами представляет собой дробно-рациональную функцию переменной p. Чтобы по- лучить функцию передачи дискретного фильтра, необходимо перейти из p-области в z-область, причем дробно-рациональный характер функции должен 48 сохраниться. Поэтому замена для переменной p должна представлять собой также дробно-рациональную функцию переменной z. Чтобы частотные харак- теристики аналогового и дискретного фильтров были связаны простой зависи- мостью, искомая замена переменной должна отображать мнимую ось в p-области на единичную окружность в z-области. В этом случае частотные ха- рактеристики аналогового и дискретного фильтров будут связаны лишь транс- формацией частотной оси, и никаких искажений «по вертикали» не будет. Простейшей из функций, удовлетворяющих перечисленным требовани- ям, является билинейное преобразование (bilinear transformation): 1 д 1 2 1 1 2 1 1 z z p F T z z − − − − = = + + Частотные характеристики аналогового а ( ) K ω ɺ и дискретного д ( ) K ω ɺ фильтров, как уже было сказано, связаны друг с другом лишь трансформаци- ей частотной оси (для наглядного сопоставления формул частотная характе- ристика дискретного фильтра представлена как функция абсолютной частоты T ω = ωɶ ): д а 2 ( ) tg 2 T K K T ω ω = ɺ ɺ , а д 2 ( ) arctg 2 T K K T ω ω = ɺ ɺ На низких частотах, когда ωT << 1, тангенс примерно равен своему ар- гументу, откуда следует 2 tg 2 T T ω ≈ ω Поэтому в области низких частот частотные характеристики аналогово- го и дискретного фильтров почти совпадают. Далее, по мере ускорения роста функции тангенса, частотная характеристика дискретного фильтра все силь- нее сжимается по горизонтали (по сравнению с аналоговым прототипом) и на частоте Найквиста, равной π/T, достигает значения, которое частотная харак- теристика аналогового фильтра имела бы на бесконечной частоте (рис. 4.2). При билинейном z-преобразовании левая половина p-плоскости отобра- жается внутрь единичной окружности на z-плоскости, поэтому синтез по ус- 49 тойчивому аналоговому прототипу даст гарантированно устойчивый дис- кретный фильтр. 0 w 0 w д /2 w | ( )| K а w | ( )| K д w Рис. 4.2. Трансформация частотной оси при билинейном z-преобразовании Для получения дискретного фильтра с заданными частотами среза необ- ходимо скорректировать частоты среза аналогового прототипа, чтобы ком- пенсировать искажения частотной оси. Так, для синтеза дискретного ФНЧ с частотой среза ω 0д аналоговый фильтр-прототип должен иметь частоту среза ω 0а , связанную с ω 0д следующим образом: 0д 0а 2 tg 2 T T ω ω = |