close

Вход

Забыли?

вход по аккаунту

?

MonakovMirolyubov

код для вставкиСкачать
Министерство образования и науки российской федерации
Государственное образовательное учреждение
высшего профессионального образования
Санкт-Петербургский государственный университет
аэрокосмического приборостроения
ОСНОВЫ ЦИФРОВОЙ
ОБРАБОТКИ СИГНАЛОВ
И МАТЕМАТИЧЕСКОЕ
МОДЕЛИРОВАНИЕ РЭС
Методические указания
к выполнению лабораторных работ
Санкт-Петербург
2011
Составители: А. А. Монаков, А. М. Миролюбов
Рецензент доктор технических наук, профессор П. Н. Петров
В методических указаниях приведены индивидуальные задания
для лабораторных работ по дисциплинам «Основы цифровой обработки сигналов» и «Компьютерные средства моделирования РЭС», а также даны краткие методические указания к их выполнению.
Предназначены для студентов технических специальностей университета, выполняющих лабораторный практикум на факультете
радиотехники, электроники и связи.
Подготовлены к публикации кафедрой бортовой радиоэлектронной аппаратуры ГУАП.
Редактор В. П. Зуева
Верстальщик А. Н. Колешко
Сдано в набор 29.12.10. Подписано к печати 28.04.11. Формат 60×84 1/16.
Бумага офсетная. Усл. печ. л. 7,2. Уч.-изд. л. 7,8. Тираж 100 экз. Заказ № 181.
Редакционно-издательский центр ГУАП
190000, Санкт-Петербург, Б. Морская ул., 67
© Санкт–Петербургский государственный
университет аэрокосмического
приборостроения (ГУАП), 2011
Методические указания
Методические указания содержат лабораторный практикум по
двум дисциплинам: основы цифровой обработки сигналов (ЦОС) и математическое моделирование РЭС. Целью лабораторного практикума
является развитие практических навыков проектирования и моделирования цифровых систем обработки сигналов при работе в среде
MATLAB и проведения исследовательской работы.
Лабораторный практикум по основам ЦОС содержит работы
№1–№6. Лабораторные работы №7–№12 выполняются в ходе изучения курса по математическому моделированию РЭС. Количество работ, выполняемых студентами, и их состав определяется учебными
планами соответствующих специальностей и факультетом, на котором они обучаются. Лабораторные работы №10–№12 носят исследовательский характер. В связи с этим перед их выполнением студент должен составить и согласовать с преподавателем план исследований по
изучаемой теме.
Студенты допускаются к выполнению лабораторной работы только
при успешной сдаче зачета по предыдущей работе и при соответствующей предварительной подготовке.
Для подготовки к работе необходимо ознакомиться с основными теоретическими сведениями по теме выполняемой лабораторной работы,
предлагаемыми примерами написания программ на языке MATLAB
и провести необходимые аналитические расчеты в соответствии с вариантом индивидуального задания. Номер индивидуального задания
сохраняется за студентом на всем протяжении лабораторного практикума.
Отчет по лабораторной работе должен содержать титульный лист,
номер и само индивидуальное задание, результаты аналитических
расчетов (если таковые предусмотрены в выполняемой лабораторной
работе), результаты выполнения написанной самостоятельно программы (численные результаты и графики), листинг программы и выводы
по работе.
В процессе подготовки к зачету следует обратить внимание на контрольные вопросы, приводимые в конце описания каждой лабораторной работы, и пользоваться литературой, перечисленной в указателе.
3
Лабораторная работа № 1
Дискретизация аналоговых сигналов
Цель работы: изучение основных свойств спектральной функции дискретных сигналов и метода ее расчета на ЦВМ.
Общие теоретические положения
Дискретным сигналом называется числовая последовательность
{s[n]}, каждому элементу которой может быть присвоен номер n,
принадлежащий множеству целых чисел Z = {…,–1,0,1,2,…}. Члены
числовой последовательности {s[n]} называются отсчетами дискретного сигнала, а множество их номеров {n} – дискретным временем.
При реализации устройств цифровой обработки сигналов дискретные сигналы получаются из аналоговых путем взятия отсчетов в равноотстоящие друг от друга моменты времени
S[n] = s(nT), n = …,–0,1,2,...,
где s(t) – аналоговый сигнал; T – интервал между соседними моментами времени. Этот процесс называется дискретизацией, а интервал T – периодом дискретизации. Взаимно однозначное превращение аналогового сигнала в дискретный, согласно теореме
Котельникова-Шеннона, возможно только в том случае, когда аналоговый сигнал имеет финитный спектр (т.е. спектральная функция отлична от нуля в конечной полосе частот). При этом возможно
обратное восстановление аналогового сигнала s(t) по соответствующему дискретному сигналу s{[n]}
s(t) = å s[n ]sinc(Fâ (t - nT )), (1.1)
n
где Fв – верхняя граница спектра аналогового сигнала,
sinc(x) = sin(πx)/(πx) – функция отсчетов.
Для выполнения (1.1) необходимо, чтобы период дискретизации
T удовлетворял критерию Котельникова-Найквиста
T£
1
.
2Fâ
(1.2)
Несоблюдение (1.2) при дискретизации приводит к так называемому наложению спектров. Рассмотрим это явление подробнее.
Для этого установим связь между спектральными функциями аналогового сигнала
4
¥
s(t)e-iω t dt
ò
Sa (ω ) =
-¥
и дискретного сигнала
Sä (Ω) =
¥
å
s[n ]e-inΩ .
n=-¥
Для этого выразим s(t) и s[n]
¥
1
2π
s(t) =
Sà (ω )eiωt dω, (1.3)
1
s[n ] =
Sä (Ω)eiΩn dΩ . 2π ò
(1.4)
ò
-¥
π
-π
Учитывая, что s[n] = s(nT), из (1.3) получим
s[n ] =
1
2π
¥
ò
Sà (ω )einωT dω. (1.5)
-¥
Сделаем в (1.5) замену переменной Ω = ωT и разобьем область интегрирования на полосы шириной 2π каждая
s[n ] =
=
¥
1
å 2πT
k=-¥
1
2πT
π
π(2k+1)
ò
π(2k-1)
æΩö
Sà çç ÷÷÷einΩ dΩ =
çè T ø
¥
æ Ω - 2kπ ö÷ inΩ
Sà çç
e dΩ.
çè T ÷÷ø
k=-¥
ò å
-π
(1.6)
Сравнивая (1.6) и (1.4), окончательно получим
Sä (Ω) =
æ Ω - 2kπ ö÷
1 ¥
÷. å Sà çç
T k=-¥ çè T ÷ø
(1.7)
Уравнение (1.7) устанавливает связь спектральных функций.
Согласно (1.7) спектр дискретизированного сигнала является суммой спектральных функций аналогового сигнала, сдвинутых друг
относительно друга на частоты, кратные частоте дискретизации
Ωд = 2π/T. Таким образом, спектральная функция сигнала s[n] является периодической с периодом Ωд функцией.
5
TU
U
4W
W
Рис. 1.1. Аналоговый сигнал и его спектр
На рис. 1.1 в качестве примера приведены графики радиоимпульса с трапециевидной спектральной функцией и сама эта функция.
На рис. 1.2 представлен график дискретизированного во времени
радиоимпульса и его спектр при несоблюдении критерия (1.2), когда частота дискретизации Ωд и верхняя частота спектра Ωв = 2πFB
удовлетворяют неравенству Ωв > Ωд–Ωв. Из рисунка видно, что в областях частот [–Ωв,–Ωд + Ωв] и [Ωв,Ωд–Ωв] верхняя боковая полоса
одного периода накладывается на нижнюю боковую полосу другого
периода функции Sa(ω). Происходит явление наложения, вследствие
6
TU
TU
4W
U
4W
W#
W½ W
#
W ½ W
#
W#
W
Рис. 1.2. Дискретный сигнал и его спектр
которого Sд(ω)≠S(ω) при |ω| ≤ Ωв. Спектр дискретной последовательности отсчетов сигнала Sд(ω) искажается по сравнению со спектральной функцией непрерывного сигнала. В то же время при выполнении условия Ωв ≤ Ωд–Ωв наложения не возникает, Sд(ω) = S(ω)
при |ω| ≤ Ωв, и восстановление аналогового сигнала по отсчетам дискретного оказывается возможным. Условие Ωв ≤ Ωд–Ωв и приводит
к критерию (1.2).
Для того чтобы условие (1.2) выполнялось на практике, входной
аналоговый сигнал системы обработки до дискретизации по времени и квантования по уровню в АЦП подвергают предваритель7
ной фильтрации с целью уменьшения полосы частот, занимаемых
сигналом. Для этого используют специальный анти-алайзинговый
(anti-aliasing) фильтр. Полосу пропускания этого фильтра выбирают из компромиссных соображений: с одной стороны для получения низкой частоты дискретизации полоса пропускания фильтра
должна быть как можно уже; с другой стороны уменьшение полосы
пропускания фильтра приводит к увеличению информационных и
энергетических потерь в системе обработки. Действительно, пусть
входной аналоговый сигнал имеет спектральную функцию Sa(ω).
Если считать, что анти-алайзинговый фильтр имеет прямоугольную частотную характеристику
ìï1, ω £ ω â
H (iω ) = ïí
,
ïïî0, ïðè äð. ω
где ωв – верхняя частота пропускания фильтра. Тогда энергетические потери в фильтре составят
Eïîò = E 1
2π
ωâ
ò
2
S (ω ) dω,
-ω â
¥
1
2
S (ω ) dω – полная энергия сигнала. Зная значение
ò
2π
-¥
энергии сигнала E и задавая величину относительных энергетических потерь Π = Eпот/E, можно определить ωв из решения следующего уравнения:
где E =
ωâ
ò
-ω â
2
S (ω ) dω = 2πE(1 - Π ).
Полученное значение ωв является исходным для расчета антиалайзингового фильтра. Период (частота) дискретизации рассчитываются после этого в соответствии с критерием КотельниковаНайквиста (1.2).
8
Индивидуальные задания к лабораторной работе
Таблица 1.1
№
Сигнал s(t)
Параметры
1
ì
2t - τ
ï
ï
,0£t£τ
ï1 s(t)= í
τ
ï
ï
ï
ï
î0, τ £ t £ T
τ = 0.5T, T = 2
2
ì
ï
τ -t
ï
,0£t£τ
ï
s(t)= í τ
ï
ï
ï
ï
î0, τ £ t £ T
τ = 0.5T, T = 2
3
ì
t
ï
ï
, 0 £ t £ τf
ï
ï
τ
f
ï
ï
ï
ïτ -t
, τ - τf £ t £ τ
s(t)= ï
í
τf
ï
ï
ï
ï
1, τ f £ t £ τ - τ f
ï
ï
ï
ï
ï
î0, τ £ t £ T
τ = 0.5T, T = 2
τ f = 0,1T
4
ì
ï æç
æ 2t - τ ÷öö÷
ï
0.5ç1 + cos ççπ
÷÷, 0 £ t £ τ
ï
ï
ç
ç
è
s(t)= í è
τ ÷ø÷÷ø
ï
ï
ï0, τ £ t £ T
ï
î
τ = 0.5T, T = 2
5
ì
æ 2t - τ ö÷
ï
÷, 0 £ t £ τ
ïï0.54 + 0.46 cos çççπ
è
s(t)= í
τ ø÷
ï
ï0, τ £ t £ T
ï
ï
î
τ = 0.5T, T = 2
6
ìï
ïï0.42 + 0.5 cos æççπ 2t - τ ÷ö÷ + 0.08 cos çæç2π 2t - τ ö÷÷, 0 £ t £ τ
çè
èç
s(t)= í
τ ø÷
τ ø÷
ïï
ïïî0, τ £ t £ T
τ = 0.5T, T = 2
7
3
ì
ï
æ 2t - τ ÷ö
ï
ï å an cos ççnπ
÷, 0 £ t £ τ
ï
çè
s(t)= í n=0
τ ÷ø
ï
ï
ï
ï
î0, τ £ t £ T
τ = 0.5T, T = 2
a0 = 0.3635819
a1 = 0.4891775
a2 = 0.1365995
a3 = 0.0106411
8
ì
ï
ïï exp éê-α t - 0.5 ùú , 0 £ t £ τ
êë
úû
s(t)= ï
τ
í
ï
ï
ï
£
£
0
,
τ
t
T
ï
î
τ = 0.5T, T = 2
α = 3.5
9
Продолжение табл.1.1
№
Сигнал s(t)
Параметры
9
2
ì
é
ï
öö÷ úù
ï
ê æ æt
ï
ï exp ê-2çççα ççç - 0.5÷÷÷÷÷÷ ú , 0 £ t £ τ
ï
øø ú
s(t)= í
ê è èτ
ï
ë
û
ï
ï
0
,
£
£
τ
t
T
ï
ï
î
τ = 0.5T, T = 2
α = 2.5
10
ì
ïeiωt , 0 £ t £ τ
s(t)= ï
í
ï
ï
î0, τ £ t £ T
τ = 0.5T, T = 2
ω = 2π * 10
11
ì
ïæ
2t - τ ÷ö iωt
ïç
÷e , 0 £ t £ τ
1ïç
ï
s(t)= íçè
τ ÷÷ø
ï
ï
ï
ï
î0, τ £ t £ T
τ = 0.5T, T = 2
ω = 2π * 10
12
ì æ
ï
æ 2t - τ ö÷ö÷ iωt
ï
÷÷ e , 0 £ t £ τ
ï0.5çç1 + cos çççπ
è
s(t)= ïí çè
τ ÷ø÷÷ø
ï
ï
ï
ï
î0, τ £ t £ T
τ = 0.5T, T = 2
ω = 2π * 10
13
2
ì
é æ æ
ù
öö÷
ïïï
t
ê
ú
ïï exp ê-2çççα ççç - 0.5÷÷÷÷÷÷ - iω tú , 0 £ t £ τ
øø
s(t)= í
è èτ
ê
ú
ï
ë
û
ï
ï
0
,
£
£
τ
t
T
ï
ï
î
τ = 0.5T, T = 2
ω = 2π * 10
14
ì
-1, 0 £ t £ 0.5τ
ï
ï
ï
ï
s(t)= í1, 0.5τ £ t £ τ
ï
ï
ï0, τ £ t £ T
ï
î
τ = 0.5T, T = 2
15
2
ìï
é æ æ
ö
öö÷ ùú
ïï 4α2 æç t
t
ê
ïï
çç - 0.5÷÷÷ exp ê-2çççα ççç - 0.5÷÷÷÷÷÷ ú , 0 £ t £ τ
ø
øø ú
s(t)= í τ è τ
ê è èτ
ïï
ë
û
ïï0, τ £ t £ T
ïî
τ = 0.5T, T = 2
α = 2.5
16
2
ìï
2
é æ æ
ö2 úù
öö÷ ùú
ï æ 2α ö êé
t
ê
2æt
ïï ççç ÷÷÷ ê1 - 4α ççç - 0.5÷÷÷ ú exp ê-2çççα ççç - 0.5÷÷÷÷÷÷ ú , 0 £ t £ τ
ï
èτ
ø ú
øø ú
s(t)= í è τ ø ê
ê è èτ
ë
û
ïï
ë
û
ïï0, τ £ t £ T
ïî
τ = 0.5T, T = 2
α = 2.5
17
2
ìï
é
é
ö÷ö úù
æt
öù
ïï
ê æ æt
ï exp ê-2çççα ççç - 0.5÷÷÷÷÷÷ ú cos êê10α ççç - 0.5÷÷÷úú , 0 £ t £ τ
ï
øø ú
èτ
øû
s(t)= í
è èτ
ê
ë
ïï
ë
û
ïï0, τ £ t £ T
ïî
τ = 0.5T, T = 2
α = 2.5
10
Окончание табл.1.1
№
Сигнал s(t)
Параметры
18
ì
1
ï
ï
,0£t£τ
ï
ï
ö2
t
ï
2æ
s(t)= í1 + α çç - 0.5÷÷÷
çè τ
ï
ø
ï
ï
ï
£
£T
0
,
τ
t
ï
ï
î
τ = 0.5T, T = 2
α =2 3
19
ì
ï æç
2ö
ï
ï
ççβ 1 -æç1 - 2t ÷ö÷ ÷÷÷
I
ï
÷
ç
0
ï
èç
τ ø÷ ÷÷
ï çççè
ø
s(t)= ï
í
,0£t£τ
ï
ï
I0 (β)
ï
ï
ï
ï0, τ £ t £ T
ï
î
τ = 0.5T, T = 2
β=6 ,
I0 ()
×
20
ì
é
ïï
æt
öù
sin ê2N π ççç - 0,5÷÷÷ú
ï
ï
êë
èτ
øúû
ï
ï
,0£t£τ
é
ù
æ
ö
s(t)= ï
í ê
çç t - 0,5÷÷ú
ï
2
N
π
ï ê
èç τ
ø÷úû
ï
ë
ï
ï
ï
0
,
τ
t
T
£
£
ï
î
– модифицированная
ф-ция Бесселя
τ = 0.5T, T = 2
N =3
Порядок выполнения и требования к отчету
Для выполнения лабораторной работы запустите из командной строки MATLAB файл интерфейса LW1V6_5.M (для версии
MATLAB 6.5) или LW1.M (для версии программы MATLAB 7.0
и выше). В соответствии с номером варианта необходимо создать
MATLAB функцию для расчета:
1) аналогового сигнала s[n];
2) спектра аналогового сигнала Sa(ω);
3) дискретного сигнала s[n];
4) спектра дискретного сигнала Sд(ω).
При этом функции s(t) и Sa(ω) рассчитываются аналитически или
численно, используя представление сигнала s(t) из табл. 1.1 и связь
между сигналом и его спектральной функцией. При невозможности
вычислить спектральную функцию аналогового сигнала аналитически, следует воспользоваться численным методом, суть которого
изложена в Прил. 1. Сигнал s[n] получается путем дискретизации
сигнала s(t). Период дискретизации задается в окне ввода интерфейса лабораторной работы (файл LW1.M или LW1V6_5.M). Другим параметром ввода является имя функции, которую вы создали. В окне
«СКО» после нажатия клавиши «Пуск» (LW1.M) или сразу после
11
выбора «радиокнопки» СИГНАЛ (LW1V6_5.M) появляется результат расчета нормированной среднеквадратической ошибки (СКО)
¥
ò
e=
-¥
2
Sa (ω )- Sä (ω ) dω
¥
ò
-¥
,
2
Sa (ω ) dω
которая является мерой разницы спектров аналогового и дискретного сигнала. СКО является функцией выбранного периода дискретизации T. В отсчете по лабораторной работе приведите график
зависимости e = e(T), построенный при различных значениях T и
объясните ход кривой в выводах по работе.
В качестве примера написания функции используйте файл
MYFUNCTION.M. Текст программы приведен далее. Рекомендуется создать копию этого файла под своим именем (!!!) в каталоге \
MATLAB\WORK для модификации программы в соответствии с
индивидуальным заданием.
Отчет по лабораторной работе должен включать:
1. Титульный лист.
2. Номер задания и формулу для сигнала из табл. 1.1.
3. Результаты теоретических расчетов.
4. Листинг программы и графики полученных результатов моделирования.
5. Выводы по работе.
12
Пример программы расчета сигнала и его спектра
function [s,s0,S,S0] = mysignal(t,f)
% Функция расчета:
% s – дискретного сигнала
% s0 – аналогового сигнала
% S – спектра дискретного сигнала
% S0 – спектра аналогового сигнала
% Аналоговый сигнал – прямоугольный видеоимпульс.
% Входные параметры функции:
% t – массив значений времени, t = (0:dt:T-dt)’,
% f – массив значений частоты, f = (-0.5/dt:df:0.5/dt-%°df)’,
%Глобальные переменные:
global T N Ti dt df
% T – интервал наблюдения сигнала
% N – количество отсчетов сигнала на интервале %°наблюдения
% Ti – длительность импульса
% dt – период дискретизации
% df – шаг расчета спектра по оси частот, df = 1/(N*dt)
% Количество подынтервалов, на которое разбивается %°период
дискретизации при численном расчете S0
K = 64;
% Общее количество подынтервалов при численном расчете %°S0
M = N*K;
% Длительность подынтервала
Dt = T/M;
% Вектор отсчетов времени при численном расчете S0
tnew = (0:M-1)’*Dt;
% Соответствующий шаг по частоте при численном расчете %°S0
Df = 1/(M*Dt);
% Вектор отсчетов частоты при численном расчете S0
fnew = (-0.5*M:0.5*M-1)’*Df;
% Задаем длительность импульса
Ti = 0.5*T;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Далее следует расчет спектров аналогового и
%°дискретного сигналов для сигнала типа
% «прямоугольный видеоимпульс»
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Обнуляем массив значений аналогового сигнала
s0 = zeros(M,1);
% Определение номеров отсчетов, удовлетворяющих условию %°0< = t< = Ti
Ind = find((0< = tnew)&(tnew< = Ti));
13
% В полученных значениях моментов времени полагаем %°сигнал равным ”1”
s0(Ind,1) = ones(length(Ind),1);
% Вычисление дискретного сигнала
s = s0(1:K:end);
% Вычисление спектра дискретного сигнала
S = fftshift(fft(s))*dt;
% Ключ: при “0” расчет S0 выполняется по аналитической %°формуле, при
“1” – численным %методом
key = 1;
switch key
case 0
% Вычисление спектра аналогового сигнала по %°аналитической формуле
S0 = Ti*sinc(fnew*Ti).*exp(-i*pi*fnew*Ti);
case 1
% Вычисление спектра аналогового сигнала
%°численным методом
S0 = Dt*fftshift(fft(s0));
end
% В массиве s0 оставляем только отсчеты, %°соответствующие вектору t
s0 = s;
%В массиве S0 оставляем только отсчеты, соответствующие %°вектору f
S0 = S0(find((-0.5/dt< = fnew)&(fnew< = (0.5-1/N)/dt)));
Контрольные вопросы
1. Что такое «дискретный сигнал» и как он описывается во времени?
2. Каковы отличительные свойства спектральной функции дискретного сигнала?
3. Каков смысл теоремы Котельникова? Для каких сигналов она
справедлива и почему?
4. Что такое явление наложения? Почему данное явление носит
негативный характер?
5. Какому условию должна соответствовать частота (период) дискретизации, чтобы исключить явление наложение?
6. Для чего нужен анти-алайзинговый фильтр?
7. Какова частотная характеристика анти-алайзингового фильтра? Как рассчитать полосу пропускания фильтра?
Литература: [1, с. 18–88; 2, c. 15–66; 3, с. 11–53, с. 94–113; 4, с. 17–86, с.
127–188; 5, с. 20–48, с. 143–185; 6, с. 7–26; 7, с. ; 8, с. 219–256]
14
Лабораторная работа № 2
Исследование линейных дискретных систем
Цель работы: Изучение методов анализа линейных дискретных
систем (цифровых фильтров) во временной и частотных областях.
Общие теоретические положения
Рассмотрим систему цифровой обработки сигналов (СЦОС) (см.
рис. 2.1), выполняющую преобразоваY<O>
Z<O>
ние входного дискретного сигнала x[n]
ª§ª
(воздействия) в выходной дискретный
сигнал y[n] (реакцию).
Рис. 2.1. Система цифровой
На рисунке x[n] = x(nT) – воздейобработки сигналов
ствие; y[n] = y(nT) – реакция; Т – период дискретизации.
В общем случае связь реакции системы и воздействия описывается операторным уравнением:
y = F{x},
где F– некоторый оператор.
Систему обработки называют линейной, если она обладает свойствами аддитивности и однородности.
Аддитивность. Система называется аддитивной, если реакция на сумму воздействий равна сумме реакций на каждое воздействие:
F{x1 + x2} = F{x1} + F{x2}.
Однородность. Система называется однородной, если умножение воздействия на весовой коэффициент соответствует реакции,
умноженной на тот же коэффициент:
F{αx} = αF{x}. Систему называют стационарной, если она обладает свойствами
инвариантности во времени:
y[n + l] = F{x[n + l]} для любого целого l. В соответствии с принципом стационарности задержка воздействия приводит к задержке реакции на то же время.
Мы будем рассматривать стационарные линейные дискретные
системы (ЛДС). На практике такие системы называются цифровыми фильтрами.
15
Импульсная характеристика. Импульсной характеристикой (ИХ) h[n]
ЛДС называется её реакция на единичный импульс u0[n] при нулевых начальРис. 2.2 Импульсная
ных условиях (см. рис. 2.2).
характеристика ЛДС
ì
ï1 , n = 0
где u0 [n ] = ïí
– единичный цифроï
ï
î0, n ¹ 0
вой импульс (цифровая дельта-функция). Нулевые начальные условия означают отсутствие реакции при отсутствии воздействия, что
отвечает принципу причинности, в соответствии с которым реакция не может возникнуть раньше воздействия: из x[n] = 0 при n<n0
следует, что y[n] = 0 при n<n0, где n0 – начальный момент (время
возникновения воздействия).
Переходная характеристика. Переходной характеристи
кой (ПХ) g[n] ЛДС называется её реакV <O>
H <O>
¤ª ция на единичный цифровой скачок
ì
ï1 , n ³ 0
при нулевых начальных
u1 [n ] = ï
Рис. 2.3. Переходная
í
ï
ï
î0, n < 0
характеристика ЛДС
условиях (см. рис. 2.3).
Описание ЛДС во временной области, формула свертки. Рассмотрим реакцию ЛДС на произвольное
Y <O>
Z <O>
воздействие x[n] (см. рис. 2.4).
¤ª Поскольку любой дискретный сигРис. 2.4. Линейная
нал x[n] может быть представлен в виде
V <O>
¤ª I<O>
дискретная система
¥
å
x[n ] =
m=-¥
u0 [n - m ]x[m ] , в силу линей-
ности системы мы получаем следующее выражение:
¥
ìï ¥
üï
F ïí å u0 [n - m ]x[m ]ïý = å F {u0 [n - m ]x[m ]}.
ïï
ïï
îm=-¥
þ m=-¥
В соответствии со свойством однородности:
¥
å
m=-¥
F {u0 [n - m ]x[m ]}=
+¥
å
m=-¥
x[m ]F {u0 [n - m ]}.
Наконец по свойству стационарности и определению импульсной
характеристики:
+¥
16
å
m=-¥
x[m ]F {u0 [n - m ]}=
+¥
å
m=-¥
x[m ]h[n - m ].
Сделав в последней сумме замену переменной l = n – m, получим
тождественную формулу
+¥
å
m=-¥
x[m ]F {u0 [n - m ]}=
+¥
å
h[m ]x[n - m ],
m=-¥
Таким образом, мы показали, что реакция ЛДС равна свертке
воздействия и ИХ:
y[n] = x[n] * h[n],
где символом « * » обозначена операция свертки:
x[n ]* h[n ] =
+¥
å
x[m ]h[n - m ].
m=-¥
Разностное уравнение. Соотношение между воздействием x[n]
и соответствующей реакцией y[n] для линейной дискретной стационарной системы в общем случае имеет вид разностного уравнения
y[n ] =
M
å
N
bm x[n - m ]- å ak y[n - k],
m=0
k=1
где bm, m = 0,...,M и an, n = 1,...,N – коэффициенты ЛДС; M, N– константы. Порядком ЛДС называют величину Q = max{M , N}.
ЛДС называется рекурсивной, если хотя бы один из коэффициентов an, n = 1,...,N разностного уравнения отличен от нуля.
Согласно разностному уравнению реакция рекурсивной ЛДС в
любой момент времени определяется:
– текущим воздействием;
– предысторией воздействия;
– предысторией реакции.
ЛДС называется нерекурсивной, если все коэффициенты an,
n = 1,...,N разностного уравнения равны нулю: an = 0, "n > 0. Таким образом, РУ нерекурсивной ЛДС имеет следующий вид:
y[n ] =
M
å bm x[n - m].
m=0
Реакция нерекурсивной ЛДС в любой момент времени определяется:
– текущим значением воздействия;
– предысторией воздействия.
Вычислим импульсную характеристику (ИХ) h[n ] для нерекурсивной ЛДС, воспользовавшись её разностным уравнением:
17
y[n ] =
M
å bm x[n - m].
m=0
По определению ИХ имеем:
h[n ] =
M
å bmu0 [n - m]= bn .
m=0
Таким образом, h[n] = bn для любого n = 0,...,M, и ИХ нерекурсивного фильтра имеет конечную длительность. Такие ЛДС называются ЛДС с конечной импульсной характеристикой (КИХ системы). ИХ рекурсивной ЛДС имеет бесконечную длительность, поэтому такие ЛДС называются системами с бесконечной импульсной
характеристикой (БИХ системы).
Устойчивость ЛДС. ЛДС называется устойчивой, если при
ограниченном воздействии и произвольных (конечных) начальных
условиях реакция ЛДС также ограничена.
Критерий устойчивости ЛДС. ЛДС является устойчивой, тогда и только тогда, когда ряд, составленный из значений ее ИХ сходится абсолютно:
¥
å h(n) < ¥.
n=0
Следствия:
1. Нерекурсивные ЛДС всегда являются устойчивыми.
2. ИХ устойчивой рекурсивной ЛДС имеет характер затухающей
функции времени.
Описание ЛДС в Z области. Основной характеристикой ЛДС в Z
области (в комплексной плоскости переменной z) является её передаточная функция H(z), которая определяется как Z-преобразование
(изображение) импульсной характеристики:
H (z) = Z {h[n ]}=
¥
å h[n]z-n .
n=0
Было установлено, что связь воздействия и реакции ЛДС определяется как свертка:
y[n] = h[n] * x[n]. Вычисляя Z-преобразование от левой и правой частей уравнения, на основании свойств Z-преобразования получим:
Z{y[n]} = Z{h[n]*x[n]} = Z{h[n]}Z{x[n]}
18
Следовательно,
Y(z) = Z{h[n]}Z{x[n]} = H(z) X(z),
где Y(z) = Z{y[n]}, X(z) = Z{x[n]}.
Y (z)
есть отношеX (z)
ние Z-преобразования реакции к Z-преобразованию воздействия.
Связь передаточной функции и разностного уравнения. Рассмотрим разностное уравнение ЛДС и вычислим Z-преобразование
от его обеих частей:
Таким образом, передаточная функция H (z) =
Z {y(n)}=
M
å
m=0
N
bm Z {x[n - m ]}- å ak Z {y[n - k]}.
k=1
Пользуясь свойствами Z-преобразования, получим
M
N
m=0
k=1
Y (z) = X (z) å bm z-m - Y (z) å ak z-k .
Приводя подобные члены, имеем
M
Y (z) =
å bm z-m
m=0
N
X (z).
-k
1 + å ak z
k=1
Отсюда получаем связь передаточной функции (ПФ) с коэффициентами разностного уравнения
M
H (z) =
å bm z-m
m=0
N
.
-k
1 + å ak z
k=1
Таким образом, передаточная функция ЛДС является дробнорациональной функцией, зависящей только от коэффициентов разностного уравнения и не зависящей от воздействия и реакции.
Как и любая дробно-рациональная функция ПФ однозначно
определяется своими нулями и полюсами. Можно показать, что
критерием устойчивости ЛДС в Z области является расположение
всех полюсов ПФ внутри единичной окружности.
19
Описание ЛДС в частотной области. Передаточная функция
(ПФ) вычисленная на единичной окружности, называется частотной характеристикой ЛДС H(iΩ):
H (iΩ) º H(eiΩ ) = H (z)
z=eiΩ
Ω Î (-π, π].
Параметр Ω Î (-π, π], который на Z плоскости равен полярному
углу, называется безразмерной или цифровой частотой. С круговой
частотой ω он связан простым соотношением
Ω = ωT,
где T – период дискретизации.
Используя выражение для передаточной функции ЛДС, получим следующее представление частотной характеристики (ЧХ) через ИХ:
H (iΩ) º H(eiΩ ) = H (z)
z=eiΩ
=
¥
å
h[m ]z-m
m=0
z=eiΩ
=
+¥
å
-imΩ
h[m]e
.
m=0
Сумма, стоящая последней, называется преобразованием Фурье
дискретного сигнала h[m]. Рассмотрим реакцию ЛДС на воздействие цифрового гармонического сигнала единичной амплитуды и
цифровой частоты Ω: x[n] = einΩ. В соответствии с формулой свертки, получим следующее выражение:
y[n] =
¥
¥
å h[m]x[n - m] = å h[m]ei(n-m)Ω =
m=0
= einΩ
m=0
¥
å h[m]e-imΩ = H(iΩ)einΩ .
m=0
Таким образом, частотную характеристику можно представить
как комплексную амплитуду, которую приобретает цифровой гармонический сигнал единичной амплитуды при прохождении ЛДС.
В экспоненциальной форме частотная характеристика (как комплексная функция) имеет следующее представление:
H(iΩ) = A(Ω)eij(Ω),
где A(Ω) и j(Ω) – вещественные функции. Модуль частотной характеристики называется амплитудно-частотной характеристикой
(АЧХ):
A(Ω) = |H(iΩ)|,
20
а фаза частотной характеристики называется фазочастотной характеристикой (ФЧХ):
j(Ω) = arg(H(iΩ)). Следовательно, значение АЧХ на частоте Ω – это амплитуда сигнала на выходе ЛДС при воздействии на входе цифрового гармонического сигнала единичной амплитуды и частоты Ω, а значение
ФЧХ на частоте Ω – это фаза выходного сигнала (сдвиг по фазе, который приобретает входной сигнал при прохождении ЛДС).
Свойства частотной характеристики:
1. ЧХ, АЧХ, ФЧХ – непрерывные функции частоты.
2. ЧХ, АЧХ, ФЧХ – периодические функции частоте с периодом
2π.
3. Если коэффициенты ПФ вещественны, то АЧХ является четной, а ФЧХ – нечетной функцией частоты.
Фазовое время задержки τф(Ω) ЛДС равно
ϕ(Ω)
τ ô (Ω) = .
Ω Групповое время задержки τгр(Ω) ЛДС определяется скоростью
изменения фазочастотной характеристики и равно
dϕ(Ω)
τ ãð (Ω) = .
dΩ Фазовое время задержки равно временному запаздыванию выходного сигнала по отношению к входному при условии, что входной сигнал – гармонический на частоте Ω. Групповое время задержки равно временному запаздыванию огибающей выходного сигнала
по отношению к огибающей входного при условии, что на вход ЛДС
подается узкополосный радиосигнала на несущей частоте Ω.
Расчет частотной характеристики. По определению частотной характеристики и на основании связи передаточной функции с
разностным уравнением ЧХ ЛДС можно представить в следующем
виде:
M
( )
F (Ω) = H eiΩ =
å bm e-m(iΩ)
m=0
N
-k iΩ
1 + å ak e ( )
º
Q (iΩ)
.
P (iΩ)
k=1
Откуда получаются выражения для амплитудно-частотной и
фазово-частотной характеристик ЛДС:
21
Re(Q)2 + Im(Q)2
A (Ω) = H eiΩ =
,
Re(P)2 + Im(P)2
( )
æ Im(Q) ö÷
æ Im(P) ö÷
÷÷ - arctg çç
÷,
φ(Ω) = arg H(eiΩ ) = arctg çç
èç Re(Q) ø÷
èç Re(P) ÷ø÷
(
)
где введены следующие обозначения: Q(iΩ) – числитель частотной
характеристики; P(iΩ) – знаменатель частотной характеристики.
Полученные выражения используются для аналитического
представления АЧХ и ФЧХ ЛДС.
Условия отсутствия искажений на выходе ЛДС. Отсутствие
искажений на выходе ЛДС означает, что сигнал может быть задержан на некоторое время k0 и его амплитуда может быть изменена в
K раз. Во временной области эти условия означают, что реакция системы на произвольное воздействие x[n] должна иметь следующий
вид:
y[n] = Kx[n – k0]. Рассмотрим условие отсутствия искажений в Z области, для чего
вычислим Z-преобразование от приведенного выражения:
Z {y[n]} = Z { Kx[n - k0 ]} = KZ {x[n - k0 ]} = Kz-k0 Z {x[n]}.
Таким образом, отсутствие искажений на выходе ЛДС с необходимостью приводит к следующему выражению для передаточной
функции:
H(z) = Kz-k0 . Это означает, что частотная характеристика равна
H(iΩ) = Ke-ik0Ω . Соответственно, АЧХ и ФЧХ имеют следующие выражения:
A(Ω) = K, j(Ω) = –k0Ω,
т. е. АЧХ должна быть постоянна, а ФЧХ должна линейно зависеть от частоты. Линейная зависимость ФЧХ от частоты приводит к
постоянному значению группового времени задержки:
d(-k0 Ω)
dφ(Ω)
τ ãð (Ω) = == k0
dΩ
dΩ
. При синтезе цифровых частотных фильтров, как правило, ставят следующие условия:
22
ïìK, ïðè Ω Î (Ωmin , Ωmax )
A (Ω) = ïí
,
ïïî 0, ïðè Ω Ï (Ωmin , Ωmax )
φ(Ω) = -k0 Ω,
означающие передачу без искажений сигналов, сосредоточенных в
полосе частот (Ωmin , Ωmax ), и подавление сигналов, сосредоточенных вне этой полосы. Однако реализация поставленных условий
(по крайней мере, первого) возможна лишь приближенно, поскольку передаточные функции ЛДС являются дробно-рациональными
и, следовательно, не могут иметь бесконечно много нулей в конечной полосе частот.
Индивидуальные задания к лабораторной работе
Таблица 2.1
№
Разностное уравнение
1
1
y[n ]= x [n ]- y[n -1]
2
2
1
1
1
y[n ]= x [n ]+ y[n -1]- y[n - 2]
2
2
3
3
y[n ]= -x [n ]+ 2x [n -1]- x [n - 2]
4
y[n ]=
5
6
2
2
x [n ]+ y[n -1]
3
3
y[n ]= x [n -1]y[n ]=
1
y[n - 2]
10
2
x [n ]- 0.25y[n -1]+ 0.25y[n - 2]
3
7
1
y[n ]= x [n ]+ 2x [n -1]+ 3x [n - 2]- y[n -1]
2
8
1
1
y[n ]= x [n ]+ y[n -1]- y[n - 2]
2
3
9
1
1
y[n ]= x [n -1]- y[n -1]+ y[n - 2]
2
3
10
y[n ]=
1
1
x [n ]- y[n - 2]
6
9
23
Окончание табл. 2.1
№
11
Разностное уравнение
y[n ]=
2
1
x [n -1]- y[n -1]
9
3
12
1
y[n ]= x [n ]- y[n -1]- y[n - 2]
5
13
1
1
y[n ]= - x [n ]+ y[n -1]
4
6
14
y[n ]= x [n ]+ 2x [n -1]+ x [n - 2]
15
1
y[n ]= x [n ]- y[n -1]- 0.75y[n - 2]
2
Порядок выполнения и требования к отчету
В соответствии с номером варианта необходимо аналитически
получить:
1. Аналитическое выражение передаточной функции ЛСД.
2. Аналитическое выражение ЧХ.
3. Аналитическое выражение АЧХ.
4. Аналитическое выражение ФЧХ.
Также необходимо создать MATLAB функцию для вычисления:
1. Импульсной характеристики ЛДС.
2. Переходной характеристики ЛДС.
3. Карты полюсов передаточной функции (исследовать устойчивость ЛДС).
4. АЧХ ЛДС.
5. ФЧХ ЛДС.
6. Группового времени задержки.
В качестве примера написания функции используйте файл
LABWORK2.M. Текст программы приведен ниже. Рекомендуется создать копию этого файла под своим именем (!!!) в каталоге \
MATLAB\WORK для модификации программы в соответствии с
индивидуальным заданием.
Отчет по лабораторной работе должен включать:
1. Титульный лист.
2. Номер задания и формулу для РУ из табл. 2.1.
3. Результаты теоретических расчетов.
4. Листинг программы и графики полученных результатов.
5. Выводы по работе.
24
Пример программы расчета характеристик ЛДС
% ОСНОВЫ ЦИФРОВОЙ ОБРАБОТКИ СИГНАЛОВ
% Листинг программы к лабораторной работе #2
%
%
%
%
%
%
Пусть РУ имеет вид
y[n] = 0.5*x[n]-0.1*x[n-1]-0.3*y[n-1]
Т.о. коэф. РУ имеют следующие значения
b0 = 0.5, b1 = -0.1, a0 = 1, a1 = 0.3
В среде MatLab коэф. a0 всегда должен быть указан
явно (a0 = 1)
clear all %Чистим рабочую область
close all %Закрываем все окна
b = [0.5 -0.1]; %Коэффициенты нерекурсивной части ЛДС
a = [1 0.3]; %Коэффициенты рекурсивной части ЛДС
N = 20; %Количество отсчетов ИХ и ПХ
n = (0:N-1)’; %Дискретное время
%1
%Единичный импульс u0
u0 = zeros(N,1);
u0(1) = 1;
%Расчет импульсной характеристики
h = filter(b, a, u0);
figure(1);
subplot(211)
stem(n,u0,’r’)
grid
xlabel(‘n’)
ylabel(‘u_0[n]’)
subplot(212)
stem(n,h,’b’)
grid
xlabel(‘n’)
ylabel(’h[n]’)
%2
%Единичный скачок u1
u1 = ones(N,1);
%Расчет переходной характеристики
g = filter(b, a, u1);
figure(2);
subplot(211)
stem(n,u1,’r’)
grid
xlabel(‘n’)
ylabel(‘u_1[n]’)
subplot(212)
stem(n,g,’b’)
25
grid
xlabel(‘n’)
ylabel(‘y[n]’)
%3
% Расчет ПФ должен быть сделан аналитически, в нашем
% примере:
% H(z) = (0.5-0.1*z^(-1))/(1 + 0.3*z^(-1))
%4
% Карта нулей и полюсов
figure(3);
zplane(b, a)
title(‘Диаграмма нулей о полюсов ЛДС’)
%5
% Расчет частотной характеристики
[H,w] = freqz(b, a, [],’whole’);
%6
% АЧХ
A = abs(H);
figure(4);
plot(w, A)
grid
xlabel(‘\Omega’)
ylabel(‘A(\Omega)’)
title(‘АЧХ’)
%7
% ФЧХ
Phi = angle(H);
figure(5);
plot(w, Phi)
grid
xlabel(‘\Omega’)
ylabel(‘\phi(\Omega)’)
title(‘ФЧХ’)
%8
% ГВЗ
[TauGR, w] = grpdelay(b, a, 512, ‘whole’);
figure(6);
plot(w, TauGR)
grid
xlabel(‘\Omega’)
ylabel(‘\tau_{gr}(\Omega)’)
title(‘Групповое время задержки’)
26
Контрольные вопросы
1. Чем отличается линейная дискретная система от нелинейной?
2. Что такое импульсная характеристика ЛДС и как ее можно
рассчитать?
3. Что такое переходная характеристика ЛДС и как ее можно
рассчитать?
4. Как связаны между собой импульсная и переходная характеристики?
5. Как можно определить коэффициент передачи ЛДС по коэффициентам разностного уравнения?
6. Каков физический смысл частотной, амплитудно-частотной и
фазочастотной характеристик ЛДС?
7. Какому условию должны удовлетворять полюсы коэффициента передачи устойчивой и реализуемой ЛДС?
8. Какой смысл имеют фазовое и групповое время задержки
ЛДС?
9. В каком случае фазовое и групповое время задержки ЛДС равны?
Литература: [1, с. 18–88; 2, c. 39–66; 3, с. 41–93, с. 114–174; 4, с. 189–
248; 5, с. 20–109; 6, с. 27–40; 8, с. 257–366]
27
Лабораторная работа №3
Автоматизированное проектирование
цифровых фильтров в среде MATLAB
Цель работы: Изучение набора программных средств SPTool пакета MATLAB для автоматизированного синтеза цифровых фильтров
Цифровые фильтры
Цифровой фильтр – это линейная дискретная система, предназначенная для извлечения цифрового сигнала из действующей на
входе фильтра смеси сигнала и помехи (см. рис. 3.1).
Y<O>V<O>W<O>
¯ÁÍÉÇ»ÇÂ
ÍÁÄÕËÉ Z<O> V<O>
Рис. 3.1. Цифровой фильтр
На рисунке введены следующие обозначения: u[n] – сигнал; v[n]
– помеха; x[n] – воздействие; y[n] – реакция цифрового фильтра.
Сигнал на выходе цифрового фильтра должен как можно точнее
соответствовать исходному сигналу u[n] и по возможности полностью подавить помеху v[n]. Поскольку степень соответствия входного u[n] и выходного y[n] сигналов может быть определена по-разному,
параметры фильтра (тип фильтра, его порядок, коэффициенты, разрядность, точность алгебраических вычислений) рассчитываются
различными методами.
Наиболее простым типом цифровых фильтров являются
частотно-избирательные фильтры, задачей которых является селекция цифровых сигналов по частоте. К таким фильтрам относятся фильтры нижних частот (ФНЧ), фильтры верхних частот (ФВЧ),
полосовой и режекторный фильтры (ПФ и РФ).
В данной лабораторной работе будут рассматривать линейные
стационарные (см. лаб. раб. № 2) частотно-избирательные фильтры.
Физическая реализуемость. Фильтр является физически реализуемым, если для него выполняется принцип причинности, т. е.
реакция фильтра не опережает соответствующее воздействие: из
условия x[n] = 0 при n<n0 следует, что y[n] = 0 при n<n0.
Как и любые линейные дискретные системы, цифровые фильтры
делятся на два класса: фильтры с конечной импульсной характери28
стикой (КИХ фильтры) и фильтры с бесконечной импульсной характеристикой (БИХ фильтры).
Проектирование цифрового фильтра – это процесс, результатом
которого является определение коэффициентов фильтра и способа
его практической реализации, параметры которого отвечают заданным техническим характеристикам. Процесс проектирования
состоит из следующих этапов:
1) задание требований к техническим характеристикам разрабатываемого устройства;
2) задание требований к параметрам фильтра;
3) синтез соответствующей передаточной функции или разностного уравнения;
4) проверка синтезированного фильтра методом математического моделирования;
5) практическая реализация на выбранной элементной базе и отладка.
На первом и втором этапах, исходя из условий применения фильтра и характеристик сигналов и помех, определяются параметры
фильтра (полосы пропускания и режекции, допустимый уровень
пульсации частотой характеристики в полосах пропускания, уровень подавления сигнала в полосах режекции). Настоящая лабораторная работа посвящена изучению третьего и четвертого этапов
процесса разработки.
Требования, предъявляемые
к частотно-избирательным цифровым фильтрам
При проектировании частотно-избирательных фильтров определяют частоту дискретизации fд и предъявляют требования к виду
амплитудно-частотной характеристики (АЧХ) в основной полосе
частот (0, fд/2). Идеальные амплитудно-частотные характеристики
частотно-избирательных фильтров представлены на рис. 3.2–3.5.
™°®
™°®
G
G½ G
Рис. 3.2. АЧХ идеального ФНЧ
G
G½ G
Рис. 3.3. АЧХ идеального ФВЧ
29
™°®
™°®
G
G½ G
G
Рис. 3.4. АЧХ идеального ПФ
G
G
G½ G
Рис. 3.5. АЧХ идеального РФ
Полосой пропускания идеального фильтра нижних частот (ФНЧ)
является область (0, f1), а полосой режекции – область (f1, fд/2).
Полосой пропускания идеального фильтра верхних частот (ФВЧ)
является область (f1, fд/2), а полосой режекции – область (0, f1).
Полосой пропускания идеального полосового фильтра (ПФ) является область (f1, f2), а полосой режекции – область (0, f1) U (f2, fд/2).
Полосой пропускания идеального режекторного фильтра (РФ)
является область (0, f1) U (f2, fд/2), а полосой режекции – область
(f1, f2).
Идеальные частотно-избирательные фильтры являются физически нереализуемыми, что приводит к необходимости пересмотра
требований к форме АЧХ физически реализуемых цифровых фильтров. Для этого вводят допустимые отклонения АЧХ синтезируемого фильтра в полосах пропускания и задерживания от АЧХ идеального фильтра. Кроме того, вводят дополнительные, так называемые
переходные полосы, расположенные между полосами пропускания
и задерживания в которых требования к АЧХ синтезируемого фильтра не предъявляют. Установленные границы допустимой АЧХ полосового фильтра, а также сама допустимая АЧХ показаны на
рис. 3.6.
™°®
D
sD
E
G G
G
G
G½ G
Рис. 3.6. Допустимая АЧХ ПФ и ее границы
30
В полосе пропускания (f2, f3) допускают отклонение АЧХ от единицы на величину ±δ. В полосах задерживания (0, f1) U (f4, fд/2) допускают отклонение АЧХ от нуля на величину ε. Дополнительно к
полосам пропускания и задерживания вводят переходные полосы
(f1, f2) U (f3, f4), в которых требования к АЧХ синтезируемого фильтра не предъявляются. Аналогичным образом задаются требования
и к другим частотно-избирательным фильтрам.
В общем случае требования к синтезируемому частотноизбирательному фильтру включают следующие параметры:
1) тип избирательности фильтра (ФНЧ, ФВЧ, ПФ и т.п.);
2) частота дискретизации fд;
3) значения граничных частот (в основной полосе частот), разделяющих полосы пропускания, задерживания, а также переходные
полосы;
4) допустимые отклонения АЧХ синтезируемого фильтра от идеальной АЧХ в полосах пропускания (δ) и режекции (ε);
5) выбор метода аппроксимации АЧХ.
Указанные требования могут предъявляться как к нормированной амплитудно-частотной характеристике A (f ) Î [0,1], так и к характеристике затухания a(f ) = 20 lg A (f ) Î (-¥,0].
Автоматизированное проектирование фильтров
В лабораторной работе рассматривается автоматизированное
проектирование цифровых фильтров на примере программных
средств Signal Processing Toolbox из состава программного комплекса MATLAB. Вызов пакета Signal Processing Toolbox (SPTool)
осуществляется с помощью команды sptool, вводимой в основном командном окне программы MATLAB, как это показано на
рис. 3.7.
В результате введенной команды на экране компьютера должно
отобразиться главное окно для работы с пакетом SPTool, представленное на рис. 3.8.
Главное окно SPTool состоит из строки меню и трех полей:
Signals, Filters и Spectra. Поля Signal и Spectra предназначены для
расчета и просмотра спектров сигналов на входе и выходе синтезированного фильтра.
Поле Filters предназначено для автоматизированного процесса
синтеза цифрового фильтра. Кнопка View предназначена для просмотра характеристик существующего фильтра, кнопка New Design
предназначена для синтеза нового фильтра, кнопка Edit Design
31
Рис. 3.7. Запуск пакета Signal Processing Toolbox
Рис. 3.8. Главное окно SPTool
предназначена для изменения характеристик существующего (выделенного) фильтра и кнопка Apply предназначена для моделирования процесса фильтрации (с помощью выделенного фильтра).
32
Рис. 3.9. Окно «Filter Design»
Рассмотрим процесс автоматизированного синтеза фильтра. Для
этого в поле Filters нажмем на кнопку New Design. При этом на
экране должно отобразиться окно «Filter Design», представленное
на рис. 3.9.
Поле «Filter» содержит имя фильтра.
Поле «Sampling Frequency» содержит значение частоты дискретизации.
Поле «Algorithm» отображает выбранный алгоритм аппроксимации:
– Equiripple FIR – оптимальный фильтрации Чебышева (КИХ);
– Least Squares FIR – метод наименьших квадратов (КИХ);
– Kaiser Window FIR – окно Кайзера (КИХ);
– Butterworth IIR – фильтр Баттерворта (БИХ);
– Chebyshev Type 1 IIR – фильтр Чебышева 1-го рода (БИХ);
– Chebyshev Type 2 IIR – фильтр Чебышева 2-го рода (БИХ);
– Elliptic IIR – фильтр Золотарева – Кауэра (эллиптический)
(БИХ);
– Pole/Zero Editor – синтез фильтра методом задания полюсов и
нулей передаточной функции в Z-области.
Флажок «Minimum Order» служит для синтеза фильтра минимального порядка, удовлетворяющего введенным требованиям.
33
Поле «Type» содержит информацию о типе избирательности
фильтра:
– lowpass – ФНЧ;
– highpass – ФВЧ;
– bandpass – ПФ;
– bandstop – РФ.
Поля «Passband» и «Stopband» служат для ввода параметров полос пропускания и задерживания соответственно. Подполя «Fp1
(Fp2)» и «Fs1 (Fs2)» служат для ввода границ полос пропускания
и режекции, а подполя «Rp» («Rs») служат для ввода максимально
(минимально) допустимого ослабления в полосе пропускания (режекции).
После ввода параметров синтезируемого фильтра необходимо нажать на флажок «Auto Design» для автоматического расчета цифрового фильтра. При успешном выполнении расчета в поле «Frequency
Response» отобразиться АЧХ синтезированного фильтра.
На рис. 3.10 представлены результаты синтеза полосового фильтра со следующими параметрами:
1. Частота дискретизации: 8000 Гц.
2. Полоса пропускания
область: (800 Гц, 1200 Гц);
максимально допустимое ослабление: 6 dB.
Рис. 3.10. Результаты синтеза полосового фильтра
34
3. Полоса пропускания
область: (0, 600 Гц) U (1400 Гц, 4000 Гц);
минимально допустимое ослабление: 40 dB.
4. Метод аппроксимации: БИХ фильтр Баттерворта.
В правой части окна в поле «Measurements» отображается порядок (Order) синтезированного фильтра. В нашем случае он оказался
равен 7.
Для анализа характеристик синтезированного фильтра необходимо закрыть окно «Filter Design» с помощью пункта меню «File\
Close» и в окне «SPTool» (см. рис. 3.8) нажать кнопку «View». При
этом на экране отобразится диалоговое окно «Filter Viewer», позволяющее просмотреть следующие характеристики синтезированного фильтра:
– АЧХ в линейном и логарифмическом масштабах (Magnitude);
– ФЧХ (Phase);
– время группового прохождения (Group Delay);
– карту нулей и полюсов (Zeros & Poles);
Рис. 3.11. Характеристики синтезированного фильтра
35
– импульсную характеристику (Impulse Response);
– переходную характеристику (Step Response).
Характеристики синтезированного фильтра представлены на
рис. 3.11.
Рис. 3.12. Окно «Export from SPTool»
Рис. 3.13. Вывод передаточной функции в командное окно MATLAB
36
Для экспортирования синтезированного фильтра в программу
MATLAB необходимо в окне «SPTool» (см. рис. 3.8) выбрать пункт
меню «File\Export…». При этом на экране отобразится окно «Export
from SPTool», представленное на рис. 3.12.
В поле «Export List» необходимо:
– выбрать имя синтезированного фильтра (в нашем примере это
«filt1»);
– установить флажок «Export Filters as TF objects»;
– нажать кнопку «Export to Workspace».
После экспортирования фильтра в программу MATLAB можно
просмотреть передаточную функцию синтезированного фильтра,
напечатав имя фильтра (filt1) в основном командном окне программы, как это показано на рис. 3.13.
После введенного имени фильтра программа отображает аналитическое выражение передаточной функции синтезированного
фильтра.
Синтез цифровых фильтров с другим типом частотной избирательности и методом аппроксимации АЧХ проводится аналогично
рассмотренному примеру синтеза полосового фильтра Баттерворта.
37
Индивидуальные задания к лабораторной работе
Таблица 3.1
№
Тип
фильтра
Параметры полосы
пропускания
Параметры полосы задерживания
Метод аппроксимации
1
ФНЧ
(0, 1000 Гц)
max ослабл.: 6 dB
(1400 Гц, 4000 Гц)
min ослабл.: 40 dB
БИХ фильтр
Баттерворта
2
ФНЧ
То же
То же
БИХ фильтр
Чебышева 1-го р.
3
ФНЧ
То же
То же
БИХ фильтр
Чебышева 2-го р.
4
ФНЧ
То же
То же
БИХ фильтр
Кауэра
5
ФВЧ
(1400 Гц, 4000 Гц)
max ослабл.: 6 dB
(0, 1000 Гц)
min ослабл.: 40 dB
БИХ фильтр
Баттерворта
6
ФВЧ
То же
То же
БИХ фильтр
Чебышева 1-го р.
7
ФВЧ
То же
То же
БИХ фильтр
Чебышева 2-го р.
8
ФВЧ
То же
То же
БИХ фильтр
Кауэра
9
ПФ
(800 Гц, 1200 Гц)
max ослабл.: 6 dB
(0, 600 Гц) U (1400
Гц, 4000 Гц)
min ослабл.: 40 dB
БИХ фильтр
Баттерворта
10
ПФ
То же
То же
БИХ фильтр
Чебышева 1-го р.
11
ПФ
То же
То же
БИХ фильтр
Чебышева 2-го р.
12
ПФ
То же
То же
БИХ фильтр
Кауэра
13
РФ
(0, 600 Гц) U
(1400 Гц, 4000 Гц)
max ослабл.: 6 dB
(800 Гц, 1200 Гц)
min ослабл.: 40 dB
БИХ фильтр
Баттерворта
14
РФ
То же
То же
БИХ фильтр
Чебышева 1-го р.
15
РФ
То же
То же
БИХ фильтр
Чебышева 2-го р.
16
РФ
То же
То же
БИХ фильтр
Кауэра
38
Таблица 3.2
№
Тип
фильтра
Параметры полосы
пропускания
Параметры полосы
задерживания
Метод аппроксимации
1
ФНЧ
(0, 1000 Гц)
max ослабл.: 6 dB
(1400м, 4000 Гц)
min ослабл.: 40 dB
Оптимальный
фильтр Чебышева
2
ФНЧ
То же
То же
Метод наименьших квадратов
3
ФНЧ
То же
То же
Окно Кайзера
4
ФНЧ
(0, 1200 Гц)
max ослабл.: 6 dB
(1400 Гц, 4000 Гц)
min ослабл.: 40 dB
Оптимальный
фильтр Чебышева
5
ФВЧ
(1400 Гц, 4000 Гц)
max ослабл.: 6 dB
(0, 1000 Гц)
min ослабл.: 40 dB
Оптимальный
фильтр Чебышева
6
ФВЧ
То же
То же
Метод наименьших квадратов
7
ФВЧ
То же
То же
Окно Кайзера
8
ФВЧ
(1400 Гц, 4000 Гц)
max ослабл.: 6 dB
(0, 1200 Гц)
min ослабл.: 40 dB
Метод наименьших квадратов
9
ПФ
(800 Гц, 1200 Гц)
max ослабл.: 6 dB
(0, 600 Гц) U (1400
Гц, 4000 Гц)
min ослабл.: 40 dB
Оптимальный
фильтр Чебышева
10
ПФ
То же
То же
Метод наименьших квадратов
11
ПФ
То же
То же
Окно Кайзера
12
ПФ
(800 Гц, 1200 Гц)
max ослабл.: 6 dB
(0, 600 Гц) U (1200
Гц, 4000 Гц)
min ослабл.: 40 dB
Окно Кайзера
13
РФ
(0, 600 Гц) U (1400
Гц, 4000 Гц)
max ослабл.: 6 dB
(800 Гц, 1200 Гц)
min ослабл.: 40 dB
Оптимальный
фильтр Чебышева
14
РФ
То же
То же
Метод наименьших квадратов
15
РФ
То же
То же
Окно Кайзера
16
РФ
(0, 600 Гц) U (1200
Гц, 4000 Гц)
max ослабл.: 6 dB
(800 Гц, 1200 Гц)
min ослабл.: 40 dB
Оптимальный
фильтр Чебышева
39
Порядок выполнения и требования к отчету
В соответствии с номером варианта необходимо синтезировать с
помощью пакета SPTool линейные КИХ и БИХ фильтры. Задание
на проектирование КИХ фильтра содержится в табл. 3.1, задание на
проектирование БИХ фильтра содержится в табл. 3.2. Во всех вариантах частота дискретизации должна быть равна 8 кГц.
Отчет по лабораторной работе должен включать:
1. Титульный лист.
2. Номер задания из табл. 3.1 и табл. 3.2.
3. Результаты синтеза:
– передаточная функция;
– график АЧХ;
– график ФЧХ;
– график времени групповой задержки;
– карта нулей и полюсов, а также заключение об устойчивости
фильтра;
– импульсная характеристика;
– переходная характеристика;
– структурная схема фильтра.
4. Выводы по работе.
Контрольные вопросы
1. Назовите основные типы частотно-избирательных фильтров.
2. Какие исходные данные необходимы для синтеза частотноизбирательного фильтра?
3. Назовите основные типы аппроксимаций АЧХ фильтров. Каковы их особенности?
4. Почему для БИХ и КИХ фильтров применяются существенно
различные методы синтеза?
5. Каковы основные характеристики цифровых фильтров?
Литература: [1, с. 89–326; 2, c. 138–199; 3, с. 114–174; 4, с. 313–372; 5,
с. 265–422; 6, с. 41–97; 8, с. 391–492]
40
Лабораторная работа № 4
Методы непараметрического спектрального анализа
случайных стационарных процессов
Цель работы: изучение классических алгоритмов непараметрического цифрового спектрального анализа случайных стационарных процессов.
Общие теоретические положения
Стационарным (в широком смысле) называется случайный
процесс ξ(t), математическое ожидание Eξ = E{ξ(t)} которого не зависит от времени, Eξ = const, а корреляционная функция (КФ)
R ξ(t1,t2) = E{(ξ(t1) – Eξ)(ξ(t2) – Eξ)} зависит только от разницы моментов времени t1 и t2, R ξ(t1,t2) = R ξ(t2 – t1).
Для стационарного случайного процесса можно математически
корректно ввести понятие спектральной плотности мощности
(СПМ) Sξ(f). Данная функция характеризует распределение средней мощности стационарного случайного процесса Pср по частотам
и ее физический смысл следует из приближенного равенства
Pср(f, f + D f) ≈ Sξ(f) Df,
(4.1)
где Pср(f, f + Df) – средняя мощность процесса, приходящаяся на полосу частот (f, f + Df). Равенство (4.1) выполняется тем точнее, чем
меньше полоса Df. Как следует из (4.1), размерность СПМ равна [Вт/
Гц]. Поэтому иногда говорят, что СПМ равна средней мощности случайного процесса в окрестности частоты f, которая приходится на
полосу Df = 1 Гц. Суммируя правую и левую части (4.1) по соседним
полосам частот, получим равенство, связывающее полную среднюю
мощность случайного процесса и СПМ,
¥
Pñð =
ò
S (f )df. (4.2)
-¥
Поскольку мощность случайного процесса – величина положительная, Sξ(f) не может принимать отрицательные значения, т.е.
для любой частоты Sξ(f)°≥°0.
Между КФ и СПМ стационарного случайного процесса ξ(t) существует взаимно однозначная связь, которую устанавливает теорема
Винера – Хинчина
¥
S (f ) =
ò
-¥
-i2πf τ
R (τ )e
¥
dτ, R (τ ) =
ò
S (f )ei2πf τ df, (4.3)
-¥
41
т.е. КФ и СПМ представляют собой пару преобразований Фурье.
Учитывая, что средняя мощность случайного процесса равна его
дисперсии, Pср = σ2, а та в свою очередь равна значению КФ в точке
τ = 0, σ2 = R ξ(0), на основании (4.2) получим, что
¥
ò
R (0) =
S (f )df, (4.4)
-¥
Данное равенство полностью соответствует второму уравнению
(4.3), если в последнее подставить τ = 0.
Равенства (4.3) позволяют рассчитать по одной из характеристик
случайного процесса (КФ или СПМ) другую. Однако на практике часто ни одна, ни другая функции не известны. Предметом спектрального анализа является оценка СПМ по наблюдаемой в ходе эксперимента на интервале [0,Tн] реализации случайного процесса ξ(t).
Оценка СПМ случайных процессов является часто встречающейся задачей при цифровой обработке. К настоящему времени
предложено большое число методов спектрального оценивания, которые можно разделить на два больших класса: непараметрические
(линейные) и параметрические (нелинейные) методы. В непараметрических методах не делается никаких предположений относительно механизма генерации наблюдаемого случайного процесса. В
параметрических методах делаются определенные предположения
о том, как наблюдаемый процесс был создан. В настоящей лабораторной работе мы ограничимся рассмотрением двух классических
методов спектрального анализа.
Метод коррелограмм. Как следует из (4.3), СПМ и КФ случайного стационарного процесса ξ(t) связаны прямым преобразованием
Фурье
¥
Sξ (f ) =
ò
Rξ (τ )e-i2πf τ dτ. (4.5)
-¥
Пусть корреляционная функция оценена в дискретные моменты
времени nT, n = 0,1,..., N (об оценивании КФ см. далее), где T – период дискретизации. Тогда оценка СПМ может быть вычислена, если
интеграл в (4.5) вычислить методом прямоугольников
42
Sˆξ (f ) = T
N-1
å
n=-N
Rˆ ξ [n ]e-i2πnfT , (4.6)
где R ξ[n] = R ξ(nT), R ξ[n] при n < 0 вычисляются на основании известного свойства корреляционной функции: R ξ[–n] = R*ξ[n], n = 1,2,….
Расчет Sξ(f) удобно производить в дискретных точках
fm = m
1
, m = 0,,2N -1. 2NT
(4.7)
Тогда
Sˆξ [m ] = T
N-1
å
n=-N
-i
Rˆ ξ [n ]e
2π
mn
2N
, m = 0,,2N -1. (4.8)
Последнее равенство целесообразно переписать в виде удобном
для использования при вычислениях алгоритма быстрого преобразования Фурье (БПФ)
Sˆξ [m ] = T
M-1
å
n=0
2π
-i mn

Rξ [n ]e M , m = 0,..., M -1, (4.9)

где M = 2N, R [n ], n = 0,, M -1 – вектор размерности 2N.
ì
ˆ
ï

ïRξ [n ], 0 £ n £ N -1
Rξ [n ] = í
.
ï
Rˆ 2N - n ], N £ n £ 2N -1
ï
ï
î ξ[
(4.10)
Отметим, что для использования алгоритма БПФ N должно быть
равно целой степени числа 2.
Пусть на спектроанализатор поступает дискретный стационарный случайный процесс ξ[n] = exp[-i(Ω0n + j)], n = 0,±1,..., где j –
случайная фаза и Ω′ = 2πf ‘T – цифровая частота. КФ такого сигнала равна R ξ[n] = exp[iΩ′n], n = 0,±1… . СПМ, соответствующая такой ФК, равна дельта-функции в точке f = f′: Sξ(f) = δ(f – f′). Однако
так будет только при бесконечном времени наблюдения процесса.
Мы же наблюдаем сигнал на конечном интервале n = 0,1,...,N–1 и,
следовательно, нам доступны в виде оценок только 2N отсчетов КФ
R ξ[n], n = –N,...,–1,0,1…,N–1. Можно сказать, что мы наблюдаем КФ
сквозь прямоугольное весовое окно шириной 2N. Это приводит к тому, что оценка СПМ перестает быть дельта-функцией. Действительно, подставляя КФ R ξ[n] в (4.8), получим
-iπ(f ¢-fm )T
Sξ [m ] = Te
sin éêë2πN (f ¢ - fm )T ùúû
sin éêë π (f ¢ - fm )T ùúû
,
(4.11)
43
где fm = m/(2N), m = 0,1,...,2N–1. Как следует из (4.11), оценка СПМ
случайного гармонического сигнала получилась в общем случае
комплексной (!). Наш результат в действительности является сверткой истинной СПМ Sξ(f) = δ(f – f′) и функции
D2N (f ) = Te-iπfT
sin[2πNfT ]
,
sin[πfT ]
(4.12)
которая называется ядром Дирихле. График модуля этой функции
от частоты f представлен на рис. 4.1. Ядро Дирихле является знакопеременной функцией, которая представляет собой один главный лепесток в окружении боковых. Максимум главного лепестка
равен 2N, а уровень боковых лепестков по отношению к максимуму
главного лепестка равен –13 дБ и не зависит от N. Ширина главного
лепестка по первым нулям равна 1/N, т.е. обратно пропорциональна количеству наблюдаемых отсчетов сигнала. Следовательно, разрешающая способность коррелограммного спектроанализатора может быть оценена как ширина окна Дирихле и равна 1/N, т.е. обратно пропорциональная длине интервала наблюдения сигнала.
Из приведенного примера также следует, что использование коррелограммного метода в изложенном виде приводит к возникновению значительных пульсаций у оценки СПМ. Причиной этих пульсаций является обрезание «хвостов» корреляционной функции при
переходе от интеграла (4.5) к конечной сумме (4.6). Поскольку это обрезание эквивалентно перемножению бесконечной по своей протяженности корреляционной функции на прямоугольное окно протя
/
%/ G
5G
Рис. 4.1. Ядро Дирихле D2N(f), N = 16
44
женности 2N, то здесь мы встречаемся с эффектом Гиббса [6, с. 76].
Для уменьшения уровня пульсаций используются
различные весо
вые окна, на которые умножаются отсчеты R [n ], n = 0,..., M -1. При
этом уровень пульсаций удается уменьшить. Однако разрешающая
способность такого спектроанализатора хуже, чем у спектроанализатора с прямоугольным весовым окном.
Метод периодограмм. Пусть в результате эксперимента получена дискретная выборка отсчетов процесса ξ[0],..., ξ[N–1]. Без потери
общности можно считать, что математическое ожидание этих отсчетов равно нулю. Вычислим для полученной выборки дискретное
преобразование Фурье (ДПФ)
Gξ [m ] =
N-1
-i
å ξ[n]e
2π
mn
N
, m = 0, N -1, (4.13)
n=0
где Gξ[m] – дискретный отсчет мгновенного спектра процесса G(f)
m
при fm =
. Вычислим TN–1|G[m]|2:
NT
2
TN-1 Gξ [m ] = TN-1
N-1 N-1
å
-i
å ξ[n]ξ* [k]e
2π
m(n-k)
N
.
(4.14)
n=0 k=0
Простой заменой индексов суммирования двойную сумму в (4.14)
можно привести к виду
N-1
2
TN-1 Gξ [m ] = T
å
n=-(N-1)
-i
Rˆ ξ [n ]e
2π
mn
N
,
(4.15)
где
1
Rˆ ξ [m ] =
N
N- m -1
å
n=0
ξ[n ]ξ éën + m ùû -
смещенная оценка КФ процесса (см. ниже). Следовательно, правая
часть равенства (4.15) является ДПФ оценки КФ и может быть принята в силу теорему Винера – Хинчина за оценку СПМ. Таким образом, оценка СПМ может быть получена с использованием периодограммы (мгновенного спектра) Gξ[m] как
2
T
Sˆξ [m ] = Gξ [m ] , m = 0,, N -1. N
(4.16)
В этом и заключается смысл метода периодограмм.
45
В теории спектрального анализа доказывается, что оценки (4.9) и
(4.16) несостоятельны, т.е. при увеличении N их дисперсия не стремится к нулю. В методе периодограмм уменьшение ошибки оценивания достигается тем, что:
– производится усреднение в скользящем окне соседних отсчетов
периодограммы (метод Даньелла);
– вся выборка разбивается на неперекрывающиеся множества
отсчетов процесса с последующим равновесным усреднением периодограмм этих множеств (метод Барлетта);
– выборка процесса разбивается на перекрывающиеся множества отсчетов процесса с последующим усреднением периодограмм
этих множеств в весовом окне (метод Уэлча).
Оценка корреляционной функции. Непараметрические методы
спектрального анализа используют оценку КФ. Пусть наблюдается случайный процесс ξ(t) и получено N его отсчетов ξ[0],..., ξ[N-1] в
дискретные моменты времени. Для оценки корреляционной функции используют два способа.
Несмещенная оценка. Оценка значения корреляционной функции в момент mT равна
Rˆ ξ [m ] =
1
N- m
N- m -1
å
ξ [n ]ξ [n + m ], (4.17)
n=0
где ξ [n ] = ξ[n ]- µˆ – центрированный отсчет процесса, т.е. отсчет, из которого вычтена оценка математического ожидания
N-1
µˆ = N-1 å ξ[n ] . Оценка (4.17) является несмещенной, т.е. ее матеn=0
матическое ожидание равно истинному значению оцениваемого параметра, E {Rˆ ξ [m ]}= Rξ [m ] . Однако эта оценка обладает одним недостатком: при увеличении m дисперсия оценки R̂ξ [m ] растет. Это
объясняется тем, что для расчета R̂ξ [m ] используется (N – m) пар
отсчетов случайного процесса. В связи с этим при увеличении m количество слагаемых в сумме (4.17) уменьшается и качество усреднения ухудшается. Поэтому на практике чаще используют другую
оценку КФ.
Смещенная оценка. В этом случае оценка вычисляется как
N- m -1
æ
mö
1
R ξ [m ] =
ξ [n ]ξ [n + m ] = ççç1 - ÷÷÷ Rˆ ξ [m ]. (4.18)
å
çè
N n=0
N ø÷
Данная оценка хотя и имеет смещение, но дисперсия ошибки
оценивания для нее меньше, чем у несмещенной оценки. Понять
46
причину этого несложно, если учесть, что смещенная оценка может
быть получена из несмещенной умножением последней на треугольное окно w[n] = (1 – |m|/N), отсчеты которого убывают с увеличением m.
Генерация случайного процесса
методом формирующего фильтра
В данной лабораторной рабо
¯ÁÍÉÇ»ÇÂ
те при выполнении потребуется
Z<O>
X<O>
ÍÁÄÕËÉ
сформировать отрезок случайного
процесса. Для этого используется
метод формирующего фильтра,
Рис. 4.2. Генерация случайного
суть которого состоит в том, что процесса методом формирующего
требуемый случайный процесс фильтра: ζ[n] – белый шум; ξ[n] –
генерируемый процесс
получается при пропускании дискретного белого гауссовского шума через линейный цифровой фильтр (см. рис. 4.2).
Несложно показать, что по окончании переходных процессов в
фильтре на выходе возникает случайный стационарный процесс,
СПМ которого равна
Sξ(Ω) = |H(iΩ)|2 Sζ(Ω),
(4.19)
где H(iΩ) – ЧХ формирующего фильтра; Sζ(Ω) = 1/T – СПМ белого
шума ζ[n]; Sξ(Ω) – СПМ генерируемого процесса; Ω = 2πfT – цифровая частота.
Белый гауссовский шум представляет собой последовательность
независимых нормально распределенных случайных величин с нулевым математическим ожиданием и дисперсией σζ2 = 1/T.
Индивидуальные задания к лабораторной работе
Таблица 4.1
№
Разностное уравнение формирующего фильтра
Тип окна
1
1
y[n ] = x[n ]- y[n -1]
2
Bartlett
2
1
1
1
y[n ] = x[n ]+ y[n -1]- y[n - 2]
2
2
3
Bartlett-Hanning
3
y[n ] = -x[n ]+ 2x[n -1]- x[n - 2]
Blackman
47
Окончание табл. 4.1
№
Разностное уравнение формирующего фильтра
Тип окна
4
2
2
y[n ] = x[n ]+ y[n -1]
3
3
Blackman-Harris
5
y[n ] = x[n -1]-
1
y[n - 2]
10
Bohman
6
2
y[n ] = x[n ]- 0.25y[n -1]+ 0.25y[n - 2]
3
Chebyshev
7
1
y[n ] = x[n ]+ 2x[n -1]+ 3x[n - 2]- y[n -1]
2
Flat Top
8
1
1
y[n ] = x[n ]+ y[n -1]- y[n - 2]
2
3
Gaussian
9
1
1
y[n ] = x[n -1]- y[n -1]+ y[n - 2]
2
3
Hamming
10
1
1
y[n ] = x[n ]- y[n - 2]
6
9
Hann
11
2
1
y[n ] = x[n -1]- y[n -1]
9
3
Kaiser
12
1
y[n ] = x[n ]- y[n -1]- y[n - 2]
5
Nuttall
13
1
1
y[n ] = - x[n ]+ y[n -1]
4
6
Parzen
14
y[n ] = x[n ]+ 2x[n -1]+ x[n - 2]
Rectangular
15
1
y[n ] = x[n ]- y[n -1]- 0.75y[n - 2]
2
Tukey
Порядок выполнения и требования к отчету
В соответствии с номером варианта необходимо по заданному в
табл. 4.1 разностному уравнению аналитически получить выражение для АЧХ формирующего фильтра. Используя пример програм48
мы оценки СПМ (см. следующий раздел), написать собственную версию, куда ввести коэффициенты формирующего фильтра, результаты расчета АЧХ и тип весового окна. Получить графики оценок КФ
и СПМ коррелограммным и периодограммным методами. Отчет по
лабораторной работе должен включать:
1. Титульный лист.
2. Номер задания и расчет АЧХ для РУ из табл. 4.1.
3. Листинг программы и графики полученных результатов.
4. Выводы по работе.
Пример программы оценки СПМ
% ОСНОВЫ ЦИФРОВОЙ ОБРАБОТКИ СИГНАЛОВ
% Листинг программы к лабораторной работе #4
%
%
%
%
%
%
Пусть РУ имеет вид
y[n] = 0.5*x[n]-0.1*x[n-1]-0.3*y[n-1]
Т.о. коэф. РУ имеют следующие значения
b0 = 0.5, b1 = -0.1, a0 = 1, a1 = 0.3
В среде MatLab коэф. a0 всегда должен быть указан
явно (a0 = 1).
clear all % Чистим рабочую область
close all % Закрываем все окна
b = [0.5 -0.1]; % Коэффициенты нерекурсивной части ЛДС
a = [1 0.3]; % Коэффициенты рекурсивной части ЛДС
win_type = ’@triang’; % Тип весового окна
M = pow2(12); % Количество отсчетов случайного процесса
n = (-M:M-1)’; % Дискретное время
f = (-M:M-1)’/(2*M); % Цифровая частота
% Inline-функция для расчета AЧХ формирующего фильтра
Aline = inline(‘abs((0.5-0.1*z1)./(1 + 0.3*z1))’,’z1’);
% Количество отчетов необходимое для завершения
% переходных процессов
M_idle = 0.125*M;
% Выборка белого гауссовского шума с единичной
% дисперсией
zeta = randn(M + M_idle,1);
% Выходной сигнал формирующего фильтра
xi = filter(b,a,zeta);
% Оценка корреляционной функции (КФ) случайного
% процесса
49
R = zeros(2*M,1); % Резервируем память под отсчеты КФ
% Вычисляем отсчеты весового окна
win = window(eval(win_type),2*M);
% Вычисление несмещенной оценки КФ
R(1:2*M-1) = xcorr(xi(M_idle + 1:end))/M;
R(2*M) = R(1); % Продолжение КФ
R = R.*win; % “Взвешиваем” оценку КФ
% Коррелограммная оценка СПМ
S_est_1 = fftshift(abs(fft(fftshift(R))));
% Мгновенный спектр случайной реализации
Xi = fftshift(fft(xi(M_idle + 1:end)));
% Периодограммная оценка СПМ
S_est_2 = real(Xi.*conj(Xi))/M;
% Расчет AЧХ формирующего фильтра
A = Aline(exp(-i*2*pi*f));
S_true = real(A.*A); % Истинная СПМ
figure(1)
plot(n,R)
grid
title(‘Оценка КФ’)
xlabel(‘n’)
ylabel(‘R[n]’)
figure(2)
plot(f,S_est_1,f,S_true)
grid
title(‘Коррелограммная оценка СПМ’)
xlabel(‘f’)
ylabel(‘S(f)’)
legend(‘Оценка СПМ’,’Истинная СПМ’)
figure(3)
plot(f(1:2:end),S_est_2,f(1:2:end),S_true(1:2:end))
grid
title(‘Периодограммная оценка СПМ’)
xlabel(‘f’)
ylabel(‘S(f)’)
legend(‘Оценка СПМ’,’Истинная СПМ’)
% Проверка энергетического баланса
[max(R) sum(S_true)/(2*M) sum(S_est_1)/(2*M)… sum(S_est_2)/M]
50
Контрольные вопросы
1. Дайте определение и назовите основные свойства спектральной плотности мощности стационарного случайного процесса.
2. Как связаны СПМ и КФ стационарного случайного процесса?
Объясните эту связь физически.
3. Определите основные операции, производимые при оценке
СПМ по методу коррелограмм?
4. Определите основные свойства ядра Дирихле.
5. Что такое эффект Гиббса? В чем причины его возникновения
при спектральном анализе? Для чего используются весовые окна?
6. Определите основные операции, производимые при оценке
СПМ по методу периодограмм?
7. Как зависит качество оценки СПМ классическими методами
спектрального анализа от размерности вектора отсчетов наблюдаемого случайного процесса?
8. Какие методы оценки КФ существуют? Назовите основные
свойства этих оценок?
9. Как в лабораторной работы моделируются отсчеты наблюдаемого случайного процесса?
Литература: [1, с. 394–483; 2, c. 381–411; 3, с. ; 4, с. 249–312; 5, с. ; 6, с.;
7, с.; 8, с. 367–390; 9, с. 71–75; 10, с. 140–213]
51
Лабораторная работа № 5
Параметрический спектральный анализ случайных
стационарных процессов: метод Юла – Уолкера
Цель работы: изучение алгоритма цифрового спектрального анализа Юла – Уолкера.
Общие теоретические положения
Предположим, что дискретный стационарный случайный процесс ξ[n] есть результат прохождения белого гауссовского шума
ζ[n] с дисперсией σ2 через цифровой формирующий рекурсивный
фильтр (см. рис. 5.1), разностное уравнение которого имеет вид
ξ[n] = ζ[n] – a1ξ[n-1] – a2ξ[n-2]…– aN ξ[n-N],
(5.1)
где an, n = 1,...,N – коэффициенты фильтра. В теории спектрального анализа случайный процесс на выходе рекурсивного цифрового
фильтра называется регрессионным процессом порядка N, а коэффициенты фильтра – коэффициентами регрессии.
Умножим левую и правую части (5.1) на ξ[m], m ≤ n и найдем математическое ожидание. Получим
N
E {ξ[n ]ξ[m ]}= E {ζ[n ]ξ[m ]}- å ak E {ξ[n - k]ξ[m ]}, (5.2)
k=1
где E{·} – оператор вычисления статистического среднего. Математические ожидания в (5.2) представляют собой авто- и взаимные
корреляции, которые равны
E{ξ[n] ξ[m]} = R ξ[n – m],
E{ζ[n] ξ[n]} = σ2,
E{ζ[n] ξ[m]} = 0, n > m. (5.3)
¯ÁÍÉÇ»ÇÂÍÁÄÕËÉ
Z<O>
X<O>
sB/
sB
sB
; s
s
;
; s
Рис. 5.1. Генерация случайного процесса методом
формирующего фильтра: ζ[n] – белый шум,
ξ[n] – генерируемый процесс
52
Тогда (5.2) можно переписать в виде
N
Rξ [n - m] + å ak Rξ [n - m - k] = σ2δmn , (5.4)
k=1
ì1, m = n
ï
– дельта Кронекера. Уравнение (5.4) можно перегде δmn = ïí
ï
ï
î0, m ¹ n
писать в следующем матричном виде:
æ
Rξ [1]
ç Rξ [0]
ççç R [-1]
Rξ [0]
çç ξ
çç


çç
çèRξ [-N ] Rξ [-N + 1]
Rξ [N ] ÷ö çæ 1 ö÷ çæσ2 ÷ö

÷ ç ÷ ç ÷÷
 Rξ [N -1]÷÷÷ çç a1 ÷÷÷ ççç 0 ÷÷
֍
÷÷ çç ÷÷÷ = çç ÷÷÷. 

÷÷ ç  ÷÷ çç 0 ÷÷
÷ç ÷ ç ÷
Rξ [0] ø÷÷ èççaN ø÷ çèç 0 ÷ø÷

(5.5)
Следовательно, если корреляционная функция R ξ(τ) случайного процесса известна, коэффициенты формирующего фильтра an,
n = 1,...,N и мощность шума σ2 могут быть найдены путем решения
системы уравнений (5.5). Уравнения системы (5.5) называются нормальными уравнениями Юла – Уолкера. Стоящая в левой части (5.5)
корреляционная матрица является теплицевой. Теплицевой называется такая матрица A, любой элемент которой удовлетворяет условию A[m,n] = A[m – n]. Поэтому элементы A[m + k,n + k] и A[m,n],
стоящие на k-й диагонали матрицы, равны: A[m + k,n + k] = A[m –
n] = A[m,n]. Таким образом, у теплицевой матрицы на всех диагоналях стоят равные элементы. Для решения системы уравнений (5.5)
с теплицевой матрицей используется алгоритм Левинсона – Дербина [10], который позволяет значительно сократить количество вычислений для нахождения обратной матрицы для матрицы положительно определенной и теплицевой.
Если коэффициенты регрессии известны, то можно рассчитать
частотную характеристику формирующего фильтра H(f). Зная H(f)
и используя известное из теории случайных процессов уравнение,
связывающее СПМ процессов на входе и выходе линейного фильтра, немедленно получаем следующее уравнение для СПМ процесса
ξ[n]
σ2
2
Sξ (f ) = H (f ) Sζ (f ) =
N
1 + å an e-i2πfTn
2
,
(5.6)
n=1
где Sζ(f) = σ2 – СПМ белого гауссовского процесса ζ[n], T – период
дискретизации.
53
Метод Юла – Уолкера. Основная идея данного метода параметрического спектрального анализа заключается в решении нормальной системы (5.5) относительно неизвестных коэффициентов
регрессии an, n = 1,...,N и мощности шума σ2 после подстановки
вместо отсчетов корреляционной функции R ξ[m], m = 0,± 1,..., ±N
их оценок, полученных по наблюдаемым отсчетам случайного процесса
1
Rˆ ξ [m ] =
N
N- m -1
å
ξ[n ]ξ[n + m ]. (5.7)
n=0
После нахождения указанным способом коэффициентов регрессии, оценка СПМ вычисляется на основании (5.6).
Метод Юла – Уолкера отличается от классических методов спектрального анализа лучшей разрешающей способностью. Причину
высокой разрешающей способности можно объяснить следующим
образом. Полученные в ходе оценивания коэффициенты регрессии
дают возможность найти отсчеты КФ процесса с номерами большими, чем N. Для этого используется уравнение (5.4)
N
Rξ [n] = - å ak Rξ [n - k], n > N. (5.8)
k=1
Таким образом, можно получить оценки значений КФ, которые
не могли быть определены при использовании классических непараметрических методов оценки. Это дает возможность оценить СПМ
процесса с большей разрешающей способностью. Действительно, если количество оцененных изложенным способом значений КФ равно M, M > 2N, то, используя коррелограммый метод оценки СПМ
(см. уравнение (4.9)), имеем
Sˆξ [m ] = T
M-1
å
n=0
2π
-i mn

Rξ [n ]e M , m = 0,..., M -1. (5.9)
Следовательно, реальная разрешающая способность построенного спектроанализатора будет равна ∆f = 1/МТ, поскольку для оценки СПМ мы используем М отсчетов КФ.
Очевидным недостатком метода Юла – Уолкера, как впрочем и
всех параметрических методов, является то, что он является в некотором смысле оптимальным только для случайных процессов,
которые получаются после прохождения белого гауссовского шума
через линейный рекурсивный (БИХ) фильтр. Не всякий случайный
процесс может быть получен таким образом. Например, гауссовский случайный процесс с коэффициентом корреляции r[m] = exp[–
54
(π∆fTm)2]; ∆f – ширина спектра, не может быть получен с помощью
формирующего фильтра, поскольку его СПМ не является дробнорациональной функцией, т.е. не может быть представлена в виде
отношения двух полиномов переменной f. Другим примером являются случайные процессы, СПМ которых в некоторых точках принимает нулевые значения (СПМ (5.6) не может принимать нулевое
значение ни при каких значениях f).
Оправданием применения метода Юла – Уолкера для спектрального анализа произвольного случайного процесса может служить
тот факт, что с увеличением порядка фильтра его коэффициенты
можно подобрать так, что СПМ (5.6) неограниченно приближается к
СПМ анализируемого процесса. Основанием для этого утверждения
является теория аппроксимации Паде, утверждающая, что любую
аналитическую функцию можно неограниченно точно аппроксимировать функцией дробно-рациональной. При этом, очевидно, встает
вопрос о выборе предполагаемого порядка формирующего фильтра
N. Исследования показывают, что если порядок выбран слишком
малым, оценка СПМ получается излишне сглаженной, если же выбран излишне большой порядок, то в оценке СПМ появляются ложные пики.
Для выбора порядка формирующего фильтра предложено несколько различных критериев [10]. В настоящее время наиболее часто используемым является информационный критерий Акаике.
Индивидуальные задания к лабораторной работе
Таблица 5.1
№
Разностное уравнение формирующего фильтра
1
1
y[n ] = x[n ]- y[n -1]
2
2
1
1
1
y[n ] = x[n ]+ y[n -1]- y[n - 2]
2
2
3
3
y[n ] = -x[n ]+ 2x[n -1]- x[n - 2]
4
2
2
y[n ] = x[n ]+ y[n -1]
3
3
5
y[n ] = x[n -1]-
1
y[n - 2]
10
55
Окончание табл. 5.1
№
Разностное уравнение формирующего фильтра
6
2
y[n ] = x[n ]- 0.25y[n -1]+ 0.25y[n - 2]
3
7
1
y[n ] = x[n ]+ 2x[n -1]+ 3x[n - 2]- y[n -1]
2
8
1
1
y[n ] = x[n ]+ y[n -1]- y[n - 2]
2
3
9
1
1
y[n ] = x[n -1]- y[n -1]+ y[n - 2]
2
3
10
1
1
y[n ] = x[n ]- y[n - 2]
6
9
11
2
1
y[n ] = x[n -1]- y[n -1]
9
3
12
1
y[n ] = x[n ]- y[n -1]- y[n - 2]
5
13
1
1
y[n ] = - x[n ]+ y[n -1]
4
6
14
y[n ] = x[n ]+ 2x[n -1]+ x[n - 2]
15
1
y[n ] = x[n ]- y[n -1]- 0.75y[n - 2]
2
Порядок выполнения и требования к отчету
В соответствии с номером варианта необходимо по заданному в
табл. 5.1 разностному уравнению создать математическую модель
формирующего фильтра и сгенерировать случайный процесс. Используя пример программы оценки СПМ (см. следующий раздел),
написать собственную версию, куда ввести коэффициенты формирующего фильтра, результаты расчета АЧХ и тип весового окна. Получить графики истинной СПМ и ее оценок при разных значениях
предполагаемого порядка формирующего фильтра. Построить график ошибки оценивания СПМ
56
∆ (N ) = max Sξ (f )- Sˆξ (f )
f
от предполагаемого порядка формирующего фильтра. Отчет по лабораторной работе должен включать:
1. Титульный лист.
2. Номер задания и РУ из табл. 4.1.
3. Листинг программы и графики полученных результатов.
4. Выводы по работе.
Пример программы оценки СПМ
% ОСНОВЫ ЦИФРОВОЙ ОБРАБОТКИ СИГНАЛОВ
% Листинг программы к лабораторной работе #5
%
%
%
%
%
%
Пусть РУ имеет вид
y[n] = 0.5*x[n]-0.1*x[n-1]-0.3*y[n-1]
Т.о. коэф. РУ имеют следующие значения
b0 = 0.5, b1 = -0.1, a0 = 1, a1 = 0.3
В среде MatLab коэф. a0 всегда должен быть указан
явно (a0 = 1).
clear all % Чистим рабочую область
close all % Закрываем все окна
b = [0.5 -0.1]; % Коэффициенты нерекурсивной
% части ЛДС
a = [1 0.3]; % Коэффициенты рекурсивной части ЛДС
ord = 1; % Порядок формирующего фильтра
M = pow2(12); % Количество отсчетов случайного
% процесса
n = (-M:M-1)’; % Дискретное время
f = n/(2*M);
% Inline-функция для расчета AЧХ формирующего
% фильтра
Aline = inline(‘abs((0.5-0.1*z1)./(1 + 0.3*z1))’,’z1’);
% Количество отчетов необходимое для завершения
% переходных процессов
M_idle = 0.125*M;
% Выборка белого гауссовского шума с единичной
% дисперсией
zeta = randn(M + M_idle,1);
% Выходной сигнал формирующего фильтра
xi = filter(b,a,zeta);
% Находим коэффициенты регрессии и мощность белого
% шума, используя алгоритм Левинсона – Дербина
% (Levinson-Durbin)
[ar,er] = aryule(xi(M_idle + 1:M),ord);
57
z = exp(i*2*pi*f);
S_est = er./abs(polyval(ar,z)).^2; % оценка СПМ
A = Aline(z); % Расчет AЧХ формирующего фильтра
S_true = real(A.*A); % Истинная СПМ
figure(1)
plot(f,S_est,’-k’,f,S_true,’:k’,’LineWidth’,2)
grid
title(‘Оценка СПМ, метод Юла-Уолкера’)
xlabel(‘f’)
ylabel(‘S(f)’)
legend(‘Оценка СПМ’,’Истинная СПМ’)
Контрольные вопросы
1. В чем различие непараметрических и параметрических методов спектрального анализа?
2. Почему для регрессионного процесса справедлива система
уравнений Юла – Уолкера?
3. Почему параметрические методы спектрального анализа имеют разрешение лучше, чем непараметрические методы?
4. Каковы недостатки параметрических методов спектрального
анализа и в чем их причина?
5. Как выбирается порядок системы уравнений Юла – Уолкера?
6. Определите основные операции, производимые при оценке
СПМ по методу Юла – Уолкера?
Литература: [1, с. 394–483; 4, с. 249–312; 8, с. 367–390; 10, с. 142–163,
с. 214–315]
58
Лабораторная работа № 6
Фильтр Калмана
Цель работы: изучение оптимального алгоритма линейной фильтрации Калмана.
Общие теоретические положения
Предположим, что дискретный стационарный случайный процесс x[n] есть результат прохождения белого гауссовского шума
v[n] с дисперсией σ2 через цифровой формирующий рекурсивный
фильтр (см. рис. 6.1), разностное уравнение которого имеет вид
x[n] = v[n] – a1x[n – 1] – a2x[n – 2]…– aN x[n – N],
(6.1)
где an, n = 1,...,N – коэффициенты фильтра. В теории спектрального анализа случайный процесс на выходе рекурсивного цифрового
фильтра называется регрессионным процессом порядка N, а коэффициенты фильтра – коэффициентами регрессии.
Представим (6.1) в матричном виде
x[n] = Ax[n – 1] + v[n],
(6.2)
где x[n] = (x[n],..., x[n – N + 1])Т – N-мерный вектор отсчетов случайного процесса на выходе фильтра; v[n] = (v[n],0,...,0)Т – отсчеты
белого гауссовского шума, действующего на входе фильтра,
æ-a1 -a2
çç
çç 1
0
çç
1
A = çç 0
çç
çç 

çç
çè 0
0
 -aN-1 -aN ö÷
÷

0
0 ÷÷÷
÷

0
0 ÷÷÷ - ÷÷


 ÷÷
÷÷

1
0 ÷ø
(6.3)
¯ÁÍÉÇ»ÇÂÍÁÄÕËÉ
Z<O>
X<O>
sB/
sB
sB
s
;
s
;
; s
Рис. 6.1. Генерация случайного процесса методом
формирующего фильтра: v[n] – белый шум,
x[n] – генерируемый процесс
59
Y<O>
I
Z <O>
X <O>
Рис.
6.2. Схема наблюдения (измерения) процесса x[n]:
w[n] – белый шум, y[n] – наблюдаемый процесс
матрица размера N×N. Отметим, что векторный процесс v[n] – вырожденный, так как все элементы его корреляционной матрицы Q
равны нулю кроме элемента Q(1,1), который равен σ2.
Допустим, что наблюдению (измерению) доступен не сам выходной процесс x[n], а процесс y[n], который образуется следующим образом (см. рис. 6.2):
y[n] = hx[n] + w[n],
(6.4)
где h – постоянный коэффициент; w[n] – белый гауссовский шум с
дисперсией R.
В теории фильтрации вектор x[n] называется вектором состояния, матрица A – матрицей перехода, белый шум v[n] – формирующим шумом, белый шум w[n] – шумом наблюдения. Уравнение (6.1)
или (6.2) называются уравнениями перехода, а уравнение (6.3) –
уравнением наблюдения.
Нашей задачей является нахождение оптимальной оценки случайного процесса x[n] по наблюдениям (измерениям) процесса y[n].
Запись уравнения (6.1) в матричном виде (6.2) позволяет использовать для оценки процесса x[n] оптимальный фильтр Калмана.
Задача оптимальной фильтрации формулируется следующим
образом. Пусть уравнения перехода и наблюдения имеют вид
x[n] = Ax[n –1] + v[n],
y[n] = Hx[n] + w[n],
(6.5)
где x[n] – вектор состояния размера N; y[n] – вектор наблюдения размера M, v[n] – формирующий белый гауссовский шум с нулевым математическим ожиданием и корреляционной матрицей Q размера
M×M; w[n] – белый гауссовский шум наблюдения с нулевым математическим ожиданием и корреляционной матрицей R; A – матрица перехода размера N×N; H – матрица наблюдения размера M×N.
Требуется по наблюдаемым отсчетам векторного процесса y[n] оптимальным образом оценить вектор состояния x[n].
Рассмотрим условную плотность распределения вероятности
(ПРВ)
60
f (xn yn ,, y1 )=
=
1
f (yn ,, y1 xn ,, x1 )f (xn ,, x1 )dxn-1 dx1. (6.6)
f (yn ,, y1 ) ò
Как видно из уравнения, наблюдений y[n] зависит только от значения вектора состояния x[n]. Поэтому
n
(
)
f (yn ,, y1 xn ,, x1 )= Õ f y j x j . j=1
(6.7)
Теперь рассмотрим уравнения перехода. Из него следует, что вектор состояния x[n] зависит только от x[n–1] и не зависит от значений вектора состояний в предшествующие моменты времени. Такие
процессы называются марковскими. Следовательно,
n-1
(
)
f (xn ,, x1 ) = f (x1 ) Õ f x j+1 x j . j=1
(6.8)
Подставляя (6.7) и (6.8) в (6.6), получим
f (xn yn ,, y1 )=
=
1
f (y n ,, y1 ) ò
n
n-1
j=1
j=1
Õ f (y j x j )f (x1 ) Õ f (x j+1 x j )dxn-1 dx1. (6.9)
Аналогично
f (xn+1 y n+1,, y1 )=
=
n+1
n
1
f y j x j f (x1 ) Õ f x j+1 x j dxn dx1. Õ
ò
f (yn+1,, y1 ) j=1
j=1
(
)
(
)
(6.10)
Сравнивая (6.9) и (6.10), последнее уравнение можно записать в
следующей рекуррентной форме:
f (xn+1 yn+1,, y1 )=
=
f (yn ,, y1 )
f (yn+1 xn+1 )f (xn+1 xn )f (xn yn ,, y1 )dxn .
f (yn+1,, y1 )ò
(6.11)
Анализируя (6.5) и (6.11), отметим следующее:
1. Дробь, стоящая перед интегралом, не зависит от переменной
xn + 1 и поэтому является постоянной при записи условной ПРВ
f(xn + 1|yn + 1,...,y1).
61
2. Уравнения (6.5) являются линейными и соответствуют преобразованию гауссовских случайных процессов v[n] и w[n], поэтому
процессы x[n] и y[n] также являются совместно гауссовскими.
3. Плотности распределения вероятности f(yn + 1|xn + 1) и
f(xn + 1|xn) являются гауссовскими ПРВ.
Пусть x – случайный гауссовский вектор с математическим ожиданием m и корреляционной матрицей K. Введем для гауссовской
ПРВ вектора x следующее обозначение:
é 1
ù
T
exp ê- (x - m) K-1 (x - m)ú . (6.12)
ëê 2
ûú
(2π)L K
Тогда с учетом сделанных замечаний рекуррентное уравнение
(6.11) можно переписать в виде
g(x; m, K) =
1
g(xn+1; mn+1, Pn+1 ) =
= C g(yn+1; Hxn+1, R) ò g(xn+1; Axn , Q)g(xn ; mn , Pn )dxn , (6.13)
где С – постоянная; mn и Pn – математическое ожидание и корреляционная матрица условной плотности f(xn |yn ,...,y1). После несложных преобразований в (6.13) получим два уравнения, которые позволяют рекуррентно вычислять математическое ожидание mn и
корреляционную матрицу Pn ,
1) mn + 1 = Amn + Pn + 1HTR–1(yn + 1 – HAmn),
2) P–1n + 1 = (APnAT + Q)–1 + HTR–1H. (6.14)
Учитывая, что оптимальной байесовской оценкой x[n + 1] является математическое ожидание апостериорной ПРВ f(xn + 1 |yn + 1
,...,y1), система (6.14) представляет собой систему рекуррентных
уравнений оценивания. Устройство, реализующее систему (6.14),
называется фильтром Калмана, по имени ученого, открывшего
данный способ оценивания (фильтрации) гауссовских случайных
процессов.
Обычно в теории калмановской фильтрации уравнения (6.14) записываются в расширенном виде, который позволяет определить
смысл входящих в (6.14) величин и упростить вычисления:
1) xn + 1|n = Axn|n ,
2) xn|n = xn|n–1 + Kn(yn + 1 – Hxn|n–1),
3) Pn + 1|n = APn|nAT + Q,
4) Kn = Pn|n–1HT(HPn|n–1HT + R)–1,
5) Pn|n = Pn|n–1 – KnHPn|n–1. (6.15)
62
Связь между обозначениями (6.14) и (6.15) следующая:
– xn|n = mn – оценка вектора x[n] по выборке y[n],...,y[1],
– xn|n–1 – экстраполированная оценка x[n] по выборке y[n1],...,y[1],
– Pn|n = Pn – корреляционная матрица ошибок оценки вектора
состояния x[n] по выборке y[n],...,y[1],
– Pn|n–1 – корреляционная матрица ошибок оценки вектора состояния x[n] по выборке y[n-1],...,y[1],
– Kn – коэффициент усиления фильтра.
Система (6.15) предпочтительнее для реализации по сравнению с
системой (6.14), поскольку в уравнение 4 системы (6.15) входит операция обращения матрицы размера M×M, в то время как уравнение
2 системы (6.14) требует обращения матрицы N×N.
Вернемся к нашей задаче. Уравнение наблюдения для рассматриваемой задачи – скалярное, поэтому M = 1 и выгода от использования системы (6.15) очевидна. Матрица наблюдения равна
H = (h,0,...,0). Шум наблюдения – белый с корреляционной матрицей R = R.
Индивидуальные задания к лабораторной работе
Таблица 6.1
№
Коэффициенты формирующего фильтра
b0 = a0
a1
a2
a3
a4
1
1
–3.5897
4.8513
–2.9241
0.66301
2
1
–3.1806
3.8612
–2.1122
0.43827
3
1
–2.7737
3.0190
–1.5048
0.28793
4
1
–2.3695
2.3140,
–1.0547
0.18738
5
1
–1.9684
1.7359
–0.72447
0.12039
6
1
–1.5704
1.2756
–0.48440
0.076197
7
1
–1.1751
0.92565
–0.31042
0.047663
8
1
–0.78210
0.67998
–0.18268
0.030119
9
1
–0.39064
0.53430
–0.84234
0.020651
10
1
0
0.48603
0
0.017665
11
1
–3.8265
5.5170
–3.5516
0.86123
12
1
–3.6076
4.9796
–3.1114
0.74187
13
1
–3.3481
4.4095
–2.6894
0.63931
14
1
–3.0523
3.8275
–2.2926
0.55124
63
Окончание табл. 6.1
№
Коэффициенты формирующего фильтра
b0 = a0
a1
15
1
16
1
17
18
19
20
a2
a3
a4
–2.7245
3.2527
–1.9256
0.47566
–2.3688
2.7029
–1.5906
0.41084
1
–1.9889
2.1949
–1.2883
0.35533
1
–1.5882
1.7444
–1.0172
0.30797
1
–1.1696
1.3665
–0.77408
0.26789
1
–0.73575
1.0762
–0.55404
0.23460
Порядок выполнения и требования к отчету
В соответствии с номером варианта необходимо по заданным в
табл. 6.1 коэффициентам разностного уравнения создать математическую модель формирующего фильтра и сгенерировать случайный процесс. Используя пример программы оптимальной фильтрации (см. следующий параграф), написать собственную версию,
куда ввести коэффициенты формирующего фильтра. Построить зависимость ошибки оценивания от отношения сигнал/шум q2 = Ps/
Pn, где Ps – мощность процесса на выходе формирующего фильтра;
Pn = R – мощность шума наблюдения. Отчет по лабораторной работе должен включать:
1. Титульный лист.
2. Номер задания и коэффициенты РУ из табл. 6.1.
3. Листинг программы и графики полученных результатов.
4. Выводы по работе.
Пример программы оценки СПМ
% ОСНОВЫ ЦИФРОВОЙ ОБРАБОТКИ СИГНАЛОВ
% Листинг программы к лабораторной работе #6
%
%
%
%
%
%
Пусть ур-ние перехода имеет вид
x[n] = 0.1*v[n] + 0.9*x[n-1]-0.2*x[n-2],
где v[n] – БГШ с дисперсией 1;
Ур-ние наблюдения имеет вид
y[n] = x[n] + c*w[n],
где w[n] – БГШ с дисперсией 1.
clear all % Чистим рабочую область
close all % Закрываем все окна
b = 1; % Коэффициенты нерекурсивной части ЛДС
a = [1 -0.9 0.2]; % Коэффициенты рекурсивной части ЛДС
64
N = 2; % Порядок формирующего фильтра
q2 = -5; % Отношение сигнал/шум
% Расчет мощности сигнала на выходе фильтра
L = 1024;
G = freqz(b,a,L,’whole’);
Ps = real(G(1:end-1)’*G(1:end-1))/L;
c = sqrt(Ps/10^(0.1*q2)); % Коэффициент усиления шума
% наблюдения
M = pow2(8); % Количество отсчетов случайного процесса
n = (-M:M-1)’; % Дискретное время
M_idle = 0.125*M; % Количество отчетов, необходимое для
% завершения переходных процессов
v = randn(M + M_idle,1); % Выборка белого гауссовского шума
% с единично дисперсией
x = filter(b,a,v); % Выходной сигнал формирующего фильтра
x = x(M_idle + 1:end); % Удаляем из процесса отрезок,
% соответствующий переходному процессу в фильтре
w = randn(M,1); % Выборка белого гауссовского шума с
% единично дисперсией
y = x + c*w; % Формируем наблюдаемый процесс
% Матрица перехода
A = [-a(2) -a(3);1 0];
% Матрица наблюдения
h = 1;
H = [h 0];
% Корреляционная матрица формирующего шума
Q = zeros(N,N);
Q(1,1) = b^2;
% Корреляционная матрица шума наблюдения
R = c^2;
% Начальные условия
x_est = zeros(size(x));
Pnn = Q;
xnn = zeros(N,1);
% Процесс фильтрации
for n = 1:M
% Time update
xnn1 = A*xnn; % Эктраполированная оценка
Pnn1 = A*Pnn*A‘ + Q; % КМ экстраполированной оценки
S = H*Pnn1*H‘ + R;
K = Pnn1*H‘*inv(S); % Коэффициент усиления
65
% Measurement update
yn = y(n);
xnn = xnn1 + K*(yn-H*xnn1); % Оценка текущего вектора
% состояния
Pnn = Pnn1-K*S*K’; % КМ оценки
x_est(n) = H*xnn; % Оценка значения процесса
end
std_est = std(x-x_est); % Ошибка оценки
n = (1:M)’;
figure(1)
plot(n,x,’-’,n,y,’:’,n,x_est,’r’)
grid
legend(‘x’,’y’,’x_{est}’)
title([‘\sigma_{e} = ’ num2str(std_est)])
figure(2)
plot(n,x-x_est,’b’)
grid
legend(‘x-x_{est}’)
title([‘\sigma_{e} = ’ num2str(std_est)])
Контрольные вопросы
1. Какая модель входного сигнала используется при синтезе
фильтра Калмана? Нарисуйте ее функциональную схему.
2. Какая модель процесса наблюдения используется при синтезе
фильтра Калмана? Нарисуйте ее функциональную схему.
3. Нарисуйте функциональную схему фильтра Калмана.
4. Почему алгоритм фильтрации Калмана по существу сводится
к двум уравнениям системы (6.14)?
5. Чему равна экстраполированная оценка?
6. Как будет работать фильтр Калмана, если по какой-либо причине отсчеты наблюдаемого процесса будут недоступны?
Литература: [11, с. 234–262; 12, с. 328–410; 13]
66
Лабораторная работа № 7
Методы генерации случайных величин
с произвольным законом распределения
Цель работы: изучение методов генерации случайных величин с
произвольным законом распределения вероятности.
Метод обратной функции
Этот метод основан на следующей теореме теории вероятностей:
если имеется случайная величина h с плотностью распределения
вероятности wh(y), то случайная величина ξ
η
ξ=
ò
wη (y)dy (7.1)
-¥
имеет равномерный закон распределения на интервале [0,1]. Действительно, найдем вероятность Pr{ ξ < x} = Fξ(x), где x – некоторое
действительное число из интервала [0,1]; Fξ(x) – интегральная функция распределения случайной величины ξ. Для этого заметим, что
интеграл, стоящий в правой части (7.1), равен интегральной функции распределения случайной величины h
y
Wη (y) =
ò
wη (y)dy (7.2)
-¥
и в силу того, что wh(y) ≥ 0, является возрастающей функцией верхнего предела y. Тогда справедлива следующая цепочка равенств:
Fξ(x) = Pr{ ξ < x} = Pr{ Wh(y) < x} = =Pr{ h < Wh–1(x)} = Wh(Wh–1(x)) = x,
(7.3)
–1
где Wh (x) – функция, обратная интегральной функции распределения Wh(x). Если x < 0, то поскольку интеграл в правой части (7.1)
не может быть отрицательным, Pr{ ξ < x} = 0. Аналогично, если x > 1, то Pr{ ξ < x} = 1, так как значение этого интеграла не может быть
больше единицы. Таким образом,
ïìï0, x < 0
ï
Fξ (x) = ïíx, 0 £ x £ 1. ïï
ïïî1, x > 1
(7.4)
Дифференцируя (7.4), получим плотность распределения случайной величины ξ
67
ïìï0, x < 0
ï
(7.5)
fξ (x) = ïí1, 0 £ x £ 1. ïï
ïïî0, x > 1
Следовательно, случайная величина ξ имеет равномерное распределение в интервале [0,1]. Это дает возможность предложить следующий алгоритм генерации случайной величины с произвольным
законом распределения:
1-й шаг. Генерируется случайная величина ξ с равномерным в
интервале [0,1] законом распределения.
2-й шаг. Искомая случайная величина h получается в результате
следующих вычислений:
h = Wh–1(ξ),
(7.6)
–1
где Wh (·) – функция, обратная интегральной функции распределения Wh(·).
Пример 1. Необходимо получить случайные числа yi с плотностью распределения вероятности wh(y) = λe–λy, y ≥ 0 и интегральной
функцией вероятности Wh(y) = 1 – e–λy, y ≥ 0.
yi
Согласно теореме xi = λ ò e-λy dy. Тогда xi = Wη (yi ) = 1 - e-λyi . На0
1
ходим обратную функцию: yi = - ln(1 - xi ). Число xi распределено
λ
равномерно на интервале [0,1]. Тогда и разность 1 – xi распределена
1
равномерно. Поэтому выражение можно упростить: yi = - ln xi .
λ
Пример 2. Необходимо получить случайные числа yi , распределенные по закону Релея. У такого случайного числа плотность распределения вероятности и интегральная функция вероятности имеют соответственно вид
2
y2
y
- 2
y - 2σ 2
wη (y) =
e
, Wη (y) = 1 - e 2σ , y ³ 0.
2
σ
Случайные числа yi можно получить путем следующего преобразования равномерно распределенных в интервале [0,1] случайных
чисел xi: yi = σ -2 ln(1 - xi ) или yi = σ -2ln(xi ) .
Недостатки рассмотренного метода заключаются в том, что
– трудно найти обратную функцию (не берется интеграл в (7.1)),
– требуется достаточный расход машинного времени.
68
Метод кусочной аппроксимации плотности
распределения вероятности (Метод Н. П. Бусленко)
Суть метода состоит в замене плотности распределения вероятности ступенчатой функцией – набором K прямоугольников, вписанных в нее и имеющих одинаковые площади (см. рис. 7.1).
Предварительно перед аппроксимацией плотность распределения вероятности подвергается усечению в хвостах на интервале
[a, b]. Площади K прямоугольников должны быть одинаковыми и
равными 1/K. Выделим прямоугольник с основанием [ak ,ak + 1], его
площадь
ak+1
ò
wη (y)dy =
ak
1
.
K
(7.7)
На основании (7.7) последовательно вычисляются значения
a = a1, a2,...,ak , ak + 1,...,aK + 1 = b, начиная с точки а и заканчивая
точкой b.
Алгоритм моделирования заключается в последовательности
следующих действий:
1. Генерируется равномерно распределенное на интервале [0,1]
случайное число ξ1.
2. С помощью этого числа определяется номер k = ](K–1) ξ1 + 1[,
где «]·[» – оператор округления до ближайшего целого. Таким образом, выделяется интервал [ak , ak + 1].
3. Генерируется следующее число ξ2, равномерно распределенное
на интервале [0,1].
4. Вычисляется случайное число h = ak + (ak + 1 – ak) ξ2. Число h
является реализацией случайной величины заданного закона распределения.
XZ
B
BkBk+
b
y
Рис. 7.1. Аппроксимация плотности вероятности
ступенчатой функцией
69
Метод Бусленко реализуется при небольших K (до 64). Достоинством является малое число операций, не зависящее от точности аппроксимации (от выбора K), так как масштабирование делается заранее до моделирования. Недостатком метода является то, что точность аппроксимации не одинакова по всей области задания функции [a,b] и зависит от величины плотности wh(y). Чем меньше wh(y)
на данном интервале, тем меньше точность, так как основание вписанного прямоугольника больше.
Метод отбора Неймана (метод отказов)
Этот метод также предполагает усечение плотности вероятности
справа и слева на некотором интервале. Случайная величина h характеризуется плотностью распределения вероятности wh(y), которая усекается на интервале [a,b].
Затем генерируются два равномерно распределенных случайных
wη (y)ли точка с коорчисла: ξ1 и ξ2 и осуществляется проверка: попадает
динатами [a + (b-a)ξ1, wh maxξ2] под кривую плотности вероятности
(см. рис. 7.2). Если это так, то запоминается первое число ξ1, которое
wη max величины h = ξ1. Критеи используется для вычисления случайной
рием отбора является очевидное неравенство
wh maxξ2 ≤ wh(ξ1). (7.8)
Пары случайных чисел, удовлетворяющие
этому условию, можξ2
но рассматривать как координаты случайных чисел на плоскости,
равномерно распределенных вдоль осей y и wh(y). Вероятность то∆y
го, что случайная точка на плоскости, попавшая под кривую wh(y),
окажется в элементарной полосе с основанием [y,y + ∆y], равна, очевидно, площади этой полосы, т.е. wh(y) ∆y. Это и есть условие, необah = a + (b – a)ξ ξимела
ходимое для того, чтобы случайная величина
1 1
заданную плотность распределения вероятности wh(y). Таким обра
XHZ
XHNBY
X
$Z
B
X
C
Рис. 7.2. К методу Неймана
70
Z
зом, алгоритм получения последовательности случайных чисел, обладающих исходной плотностью, может быть сформулирован следующим образом:
1. Из исходной совокупности равномерно распределенных на интервале [0,1] чисел выбираем пары ξ1, ξ2.
2. Для этих чисел осуществляется проверка неравенства (7.8).
3. Если неравенство (7.8) справедливо, то переходим к шагу 4. В
противном случае к шагу 1.
4. Если неравенство выполняется, то очередное число определяется согласно следующему соотношению: h = a + (b–a)ξ1.
Описанная процедура отбора случайных чисел может потребовать существенного числа вычислений, в основном за счет вычисления правой части неравенства (7.4). Кроме того, не все пары чисел ξ1,
ξ2 будут удовлетворять (7.8) и, следовательно, некоторая часть этих
пар будет отброшена. Это приводит к дополнительным затратам машинного времени.
Индивидуальные задания к лабораторной работе
Таблица 7.1
№
Распределение
1
Хи-квадрат с n
степенями свободы
2
Логарифмическое и нормальное
3
Хи-квадрат с n
степенями свободы
4
Симпсона
Плотность распределения wh(y)
1
n/2
2
Γ (n / 2)
yn/2-1e-y/2 ,0 < y < ¥
ìï ln y - m 2 üï
) ïï
ï (
exp ïíý,0 < y < ¥
2
ï
ïï
2σ
ïîï
2πσ2 y2
þï
1
1
n/2
2
Γ (n / 2)
yn/2-1e-y/2 ,0 < y < ¥
0,-¥ < y < a
ïìï
ïï
a+b
ïï 4(y - a)
,a < y <
ïï
2
2
ïï (b - a)
í
ïï 4(b - y) a + b
ïï
,
<y<b
ïï (b - a)2 2
ïï
ïïî
0, b < y < ¥
Параметры
распределения
n = 2
m = 1
σ = 0,3
n = 1
a = 0
b = 2
71
Продолжение табл. 7.1
№
Распределение
5
Хи-квадрат с n
степенями свободы
Хираспределение
с n степенями
свободы
Хираспределение
с n степенями
свободы
Гаммараспределение
6
7
8
9
Плотность распределения wh(y)
1
n/2
2
Γ (n / 2)
yn/2-1e-y/2 ,0 < y < ¥
1
n/2-1
2
Γ (n / 2)
1
n/2-1
2
Гаммараспределение
Γ (n / 2)
β
α +1
β
α +1
2
yn-1e-y /2 ,0 < y < ¥
2
yn-1e-y /2 ,0 < y < ¥
1
Γ (α + 1)
1
Γ (α + 1)
Параметры
распределения
n = 4
n = 1
n = 4
yα e-y/β ,0 < y < ¥
α = 2
β = 1
yα e-y/β ,0 < y < ¥
α = 3
β = 1
10
Накагами (m –
распределение)
m
ì m 2ï
ü
ï
2 æç m ÷ö 2m-1
exp ï
í- 2 y ï
ý,0 < y < ¥
çç 2 ÷÷÷ y
ï
ï
Γ (m) è σ ø
ï σ
ï
î
þ
m = 0,5
σ = 1
11
Накагами (m –
распределение)
m
ì m 2ï
ü
ï
2 æç m ÷ö 2m-1
exp ï
í- 2 y ï
ý,0 < y < ¥
çç 2 ÷÷÷ y
ï
ï
Γ (m) è σ ø
ï σ
ï
î
þ
m = 1
σ = 1
12
Накагами (m –
распределение)
m
ì m 2ï
ü
ï
2 æç m ÷ö 2m-1
exp ï
í- 2 y ï
ý,0 < y < ¥
çç 2 ÷÷÷ y
ï
ï
Γ (m) è σ ø
ï σ
ï
î
þ
m = 2
σ = 1
13
Коши
1
t
, y <¥
π t2 + (y - m)2
t = 1
m = 1
14
Релея
ïìï y2 ïüï
exp
í- 2 ý,0 < y < ¥
ïï 2σ ïï
σ2
î
þ
σ = 1
15
Релея-Райса
ïì y2 + m2 ïüï æ my ö÷
÷,0 < y < ¥
expïíI çç
2
2 ý
ïï
ïï 0 çè σ2 ÷÷ø
σ
2
σ
î
þ
m = 1
σ = 1
16
Максвелла
72
y
y
2
π
ìï y2 üï
ï,0 < y < ¥
exp ïí3
2ý
ï
ï
σ
îï 2σ þï
y2
σ = 1
Окончание табл. 7.1
Параметры
распределения
№
Распределение
Плотность распределения wh(y)
17
Стьюдента
(t-распределение)
с m-степенями
свободы
æ m + 1ö÷
m+1
Γ ççç
÷
è 2 ÷ø æç
y2 ö÷÷ 2
ç1 + ÷
, y <¥
æ m ö çç
m ÷ø÷
πmΓ ççç ÷÷÷ è
è2ø
m = 2
18
Стьюдента
(t-распределение)
с m-степенями
свободы
æ m + 1ö÷
m+1
Γ ççç
÷
è 2 ø÷ æç
y2 ÷ö÷ 2
ç1 + ÷
, y <¥
æ m ö çç
m ÷÷ø
πmΓ ççç ÷÷÷ è
è2ø
m = 4
19
Эрланга k-го
порядка
λ k+1 k -λy
y e
,0 < y < ¥
Γ (k + 1)
λ = 2
k = 0
20 Эрланга k-го порядка
λ k+1 k -λy
y e
,0 < y < ¥
Γ (k + 1)
λ = 2
k = 1
21 Эрланга k-го порядка
λ k+1 k -λy
y e
,0 < y < ¥
Γ (k + 1)
λ = 2
k = 3
{
}
α = 1.5
c = 1
{
}
α = 2
c = 1
cα yα-1 exp -cyα ,0 < y < ¥
{
}
α = 3
c = 1
1
exp {D cos y},-π < y < π
2π I0 (D )
D = 3
22
Вейбулла
cα yα-1 exp -cyα ,0 < y < ¥
23
Вейбулла
cα yα-1 exp -cyα ,0 < y < ¥
24
Вейбулла
25
Мизеса
Порядок выполнения и требования к отчету
В соответствии с номером варианта необходимо по заданной в
табл. 6.1 ПРВ сгенерировать выборку случайных чисел двумя (как
минимум!) из рассмотренных методов. Отказ от использования метода обратной функции должен быть обоснован.
В качестве примера используйте специально написанную в качестве примера программу (см. следующий параграф). Рекомендуется создать копию этой программы под своим именем в каталоге \
73
MATLAB\WORK для последующей модификации в соответствии с
индивидуальным заданием.
Отчет по лабораторной работе должен включать:
1. Титульный лист.
2. Номер задания, заданную ПРВ и результаты ее теоретического исследования, обоснование выбранных методов генерации случайных чисел.
3. Листинг программы и графики полученных результатов моделирования.
4. Выводы по работе.
Пример программы генерации случайных чисел
с экспоненциальным законом распределения
%
%
%
%
%
ОСНОВЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ РЭС
Листинг программы к лабораторной работе #7
Программа моделирования случайных величин с
показательным законом распределения
f(x) = lambda*exp(-lambda*x), x > 0
% Закрываем все открытые окна вывода
close all
% Очищаем рабочую область (WORKSPACE)
clear
% Создаем INLINE-объект fline, т.е. строку символов,
% описывающую последовательность операций при
% вычислении плотности распределения вероятности
% f(x) = P1*exp(-P1*x), x > 0;
% P1 = lambda – параметр з-на.
% Знак “-” в начале символьной строки необходим для
% того, чтобы можно было в дальнейшем использовать
% встроенную функцию FMIN для нахождения
% максимума f(x)
fline = inline(‘-P1*exp(-P1*x)’,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
метод обратной функции
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% P1 – параметр распределения
P1 = 1;
% N – количество случайных величин, которое необходимо
% сгенерировать
N = 500;
% xi – массив случайных чисел с равномерным в интервале
% [0,1] распределением
xi = rand(N,1);
% eta – массив случайных чисел с экспоненциальным
74
% законом распределения
% Стоящее в правой части равенства выражение – обратная
% функция для интегральной функции распределения
% F(x) = 1-exp(-lambda*x), x > 0
eta = -1/P1*log(xi);
% Три нижеследующие оператора соответствуют обработке
% полученной реализации массива случайных величин eta,
% которая заключается в оценке плотности распределения.
% M – количество интервалов разбиения области значений
% генерируемой случайной величины, которое будет
% использоваться при построении гистограммы и оценке
% плотности распределения.
M = 40;
% f_est1 – оценка плотности распределения
% x – массив чисел, соответствующих серединам
% интервалов разбиения области значений генерируемой
% случайной величины
% Вычисление гистограммы
[f_est1,x] = hist(eta,M);
% Оценка плотности распределения
f_est1 = f_est1/(N*(x(2)-x(1)));
% Вычисление модельной функции распределения
f = -feval(fline,x,P1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
метод Неймана
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a – начало области значений генерируемой случайной
% величины
a = 0;
% b – конец области значений генерируемой случайной
% величины
b = 10;
% Определение точки, где плотность распределения
% генерируемой случайной величины достигает своего
% максимально значения.
% Используется встроенная функция FMIN.
Xmax = fmin(fline,a,b,[],P1);
% Fmax – максимум плотности распределения генерируемой
% случайной величины
Fmax = -feval(fline,Xmax,P1);
% n – количество удачных попыток при реализации метода
% Неймана
n = 1;
% nn – общее количество попыток при реализации метода
% Неймана
nn = 1;
% Цикл, реализующий алгоритм генерации Неймана
while n< = N
xi1 = a + (b-a)*rand(1,1);% xi1 – 1-е случайное число
75
xi2 = Fmax*rand(1,1);
% xi2 – 2-е случайное число
% Проверка условия xi2 < = f(xi1)
if xi2< = -feval(fline,xi1,P1)
eta(n) = xi1; % Удачная попытка
n = n + 1;
% Увеличиваем число удачных попыток
end;
nn = nn + 1;
% Увеличиваем общее число попыток
end;
% Вывод на экран nn и n
[nn n]
% Обработка полученной реализации
% f_est2 – оценка плотности распределения
f_est2 = hist(eta,x);
f_est2 = f_est2/(N*(x(2)-x(1)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
метод Бусленко
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% K – количество интервалов разбиения области значений
% генерируемой случайной величины, которое будет
% использоваться при реализации алгоритма Бусленко.
K = 100;
% p – вероятность попадания генерируемой случайной
% величины в один из К интервалов
p = 1/K;
% X – массив точек, соответствующих концам интервалов
X(1) = a;
% Цикл рекуррентного вычисления элементов массива X
for k = 2:K + 1
X(k) = -1/P1*log(exp(-P1*X(k-1))-p);
end;
% Цикл генерации случайных чисел с экспоненциальным
% законом распределения
for n = 1:N
xi1 = rand(1,1);% xi1 – 1-е случайное число
% k – номер случайно выбранного интервала
k = round((K-1)*xi1) + 1;
xi2 = rand(1,1);% xi2 – 2-е случайное число
% Вычисление случайной величины
eta(n) = X(k) + xi2*(X(k + 1)-X(k));
end;
% Обработка полученной реализации
% f_est3 – оценка плотности распределения
f_est3 = hist(eta,x);
f_est3 = f_est3/(N*(x(2)-x(1)));
%
76
% Построение графиков:
% f
– заданная плотность распределения, синий
цвет
% f_est1 – оценка плотности по совокупности,
полученной методом обратной функции,
зеленый цвет
% f_est2 – оценка плотности по совокупности,
%
полученной методом Неймана, красный цвет
% f_est3 – оценка плотности по совокупности,
%
полученной методом Бусленко, голубой цвет
% открываем первое окно вывода
figure(1)
% строим графики
plot(x,f,x,f_est1,x,f_est2,x,f_est3)
% вывод названия графика
title(‘Результаты моделирования’)
% метка оси абсцисс
xlabel(‘x’)
% метка оси ординат
ylabel(‘f(x)’)
% вывод сетки
grid
% вывод легенды
legend(‘истинная прв’,…
‘мет. Обр. ф-ции’,…
‘мет. Неймана’,…
‘мет. Бусленко’)
%
%
Контрольные вопросы
1. Определите вероятностный смысл плотности распределения
вероятности произвольной случайной величины. Каковы основные
свойства этой функции?
2. Определите вероятностный смысл интегральной функции распределения вероятности произвольной случайной величины. Каковы основные свойства этой функции?
3. Вычислите математическое ожидание и дисперсию случайной
величины с равномерным в интервале [0,1] распределением.
4. Определите основные шаги алгоритма метода обратной функции. Перечислите достоинства и недостатки метода.
5. Определите основные шаги алгоритма метода Неймана. Перечислите достоинства и недостатки метода.
6. Определите основные шаги алгоритма метода Бусленко. Перечислите достоинства и недостатки метода.
Литература: [7, с. 12–54; 9, c. 14–29]
77
Лабораторные работы № 8 и 9
Моделирование гауссовских случайных процессов
с заданными спектральными свойствами
Цель работы: изучение методов генерации случайных стационарных нормальных процессов.
Метод дискретного преобразования Фурье
Постановка задачи: требуется создать отрезок реализации комплексного гауссовского случайного процесса ξ(t) длительности Тн,
если известна его корреляционная функция (КФ) Kξ(τ).
В соответствии с теоремой Винера – Хинчина спектральная плотность мощности (СПМ) случайного стационарного процесса Sξ(ω)
равна
¥
Sξ (ω ) =
ò
Kξ (τ )e-iωτ dτ. (8.1)
-¥
Предположим, что процесс ξ(t) поступает на вход узкополосного
фильтра с равномерной в полосе [ω¢ – 0.5∆ω, ω¢ + 0.5∆ω] частотной
характеристикой (рис. 8.1), т.е.
ì
ï1, ω - ω ¢ £ ∆ω 2
H (iω ) = ï
,
í
ï
0, ω - ω ¢ > ∆ω 2
ï
ï
î
где ω¢ и ∆ω – средняя частота и ширина полосы пропускания фильтра. Тогда средняя мощность сигнала на выходе фильтра будет равна
1
Pη =
2π
¥
ò
-¥
1
Sξ (ω ) H (iω ) dω =
2π
2
X U
) JW ω ¢+ ∆ω /2
ò
Sξ (ω )dω. ω ¢-∆ω /2
H U
) JW $W
W
W¢
Рис. 8.1. Фильтр с равномерной
в полосе [ω¢ – 0.5∆ω, ω¢ + 0.5∆ω]
частотной характеристикой
78
(8.2)
Если интервал ∆ω настолько мал, что функция Sξ(ω) мало меняется в полосе частот [ω¢ – 0.5∆ω, ω¢ + 0.5∆ω], то будет выполняться
примерное равенство
Pη » Sξ (ω ¢)
∆ω
,
2π
(8.3)
причем (8.3) соблюдается тем точнее, чем меньше полоса ∆ω. Последнее уравнение определяет физический смысл функции Sξ(ω):
СПМ описывает распределение средней мощности случайного процесса ξ(t) по частотам гармонических колебаний, входящих в его состав. Такая интерпретация СПМ и равенства (8.3) дает возможность
предложить следующий метод генерации процесса ξ(t).
Разобьем частотную область на соприкасающиеся между собой
полосы одинаковой ширины ∆ω (см. рис. 8.2). Средняя частота m-ой
полосы равна ωm = (m + 0,5)∆ω. Пусть величина ∆ω взята настолько малой, что (8.3) выполняется с высокой точностью. Тогда сумму гармоник случайного процесса ξ(t), попадающих в m-ю полосу,
вследствие малости ∆ω, можно заменить на одно гармоническое колебание частоты ωm, которое имеет случайные амплитуду Um и фазу ϕm. При этом для всего процесса будет справедливо следующее
представление:
ξ (t) =
¥
å
m=-¥
Um exp {i[ωm t + ϕm ]}=
¥
å
m=-¥
xm eiωmt , (8.4)
где xm = Um ei ϕm – комплексная амплитуда m-й гармоники.
Из теории вероятностей известно, что сумма произвольного числа гауссовских случайных величин также является гауссовской
случайной величиной. Следовательно, поскольку процесс ξ(t) – гауссовский, т.е. в любой момент времени t¢ случайная величина ξ(t¢)
4XW
W
N $W
N $W
Рис. 8.2. Разбиение области определения
СПМ Sξ(ω) на полосы шириной ∆ω
79
распределена по нормальному закону, гармонические составляющие в правой части (8.4) также должны быть гауссовскими случайными процессами. Для этого достаточно, чтобы случайные величины Um и ϕm при различных индексах m были независимы и имели
соответственно плотность распределения Релея (Um) и равномерную
в интервале [–π,π] плотность (ϕm). Это требование эквивалентно тому, что амплитуды xm должны быть комплексными гауссовскими
случайными величинами, дисперсии которых Dm равны средней
мощности процесса ξ(t), приходящейся на m–ю полосу:
∆ω
(8.5)
,
2π
где треугольные скобки означают усреднение по ансамблю.
Если процесс ξ(t) имеет ограниченный по частоте спектр, т.е.
Sξ(ω) = 0 при |ω| > Ωmax , где Ωmax – максимальная частота спектра,
то бесконечный ряд (8.4) можно заменить конечной суммой
Dm = xm
ξ (t) =
2
= Sξ (ωm )
M /2-1
å
m=-M /2
xm eiωmt , где M ≈ 2Ωmax/∆ω – количество гармоник, необходимое для моделирования случайного процесса. Тогда дискретные отсчеты процесса
ξ[n] = ξ(nТ) во времени, взятые в соответствии с теоремой Котельникова с периодом дискретизации Т ≤ π/Ωmax могут быть получены как
ξ[n ] =
M /2-1
å
m=-M /2
i
2π
xm e N
nm
(8.6)
где ∆ωТ = 2π/N. Таким образом, отсчеты случайного процесса представляют собой дискретное преобразование Фурье (ДПФ) последовательности {xm}. Соответствующим выбором ∆ω и Т можно сделать
N равным целой степени числа 2 (N = 2L ≥ M), где L – целое число.
Это дает возможность использовать алгоритм быстрого преобразования Фурье (БПФ) для вычисления суммы (8.7).
Окончательно алгоритм генерации случайного процесса методом
дискретного преобразования Фурье может быть представлен последовательностью следующих шагов:
1-й шаг. На основании заданной КФ (или СПМ) процесса выбираются ∆ω, Т и M так, чтобы N ( = 2π/∆ωТ ≥ M) было целой степенью
числа 2.
2-й шаг. Генерируется массив из M независимых комплексных
случайных чисел zm с нулевым математическим ожиданием и единичной дисперсией.
80
3-й шаг. Вычисляются комплексные амплитуды xm = Dm zm ,
где Dm определяются в соответствии с (8.5).
4-й шаг. Вычисляются дискретные отсчеты случайного процесса
ξ[n] на основании (8.7) с использованием алгоритма БПФ.
Метод формирующего фильтра
Генерация нормального случайного процесса методом формирующего фильтра основывается на двух положениях теории случайных процессов:
1. Результатом произвольного линейного преобразования гауссовского случайного процесса является также случайный гауссовский процесс.
2. Спектральные плотности мощности случайных процессов на
входе и выходе линейного фильтра с частотной характеристикой
H(iω) связаны соотношением
2
Sâûõ (ω ) = H (iω ) Sâõ (ω ).
Если предположить, что на вход фильтра поступает процесс типа белого шума (Sбш(ω) = N0/2, |ω|< ∞), то СПМ процесса на выходе
будет равна
Sη (ω ) =
N0
2
H (iω ) . 2
(8.7)
Следовательно, для формирования гауссовского случайного процесса с заданной СПМ Sξ(ω) необходимо на вход фильтра с частотной
характеристикой, модуль которой равен
H (iω ) = Sξ (ω ) (8.8)
подать белый шум с единичной спектральной плотностью N0/2 = 1.
При этом сигнал на выходе фильтра будет иметь нормальное распределение в силу п. 1, так как фильтрация – операция линейная.
Фильтр, амплитудно-частотная характеристика (АЧХ) которого
удовлетворяет (8.8), называется формирующим.
Равенство (8.8) дает возможность вычислить лишь АЧХ и не
определяет фазочастотную характеристику (ФЧХ) фильтра. Поэтому этот вопрос остается открытым и требует дополнительного исследования.
Допустим, что ФЧХ фильтра ϕ(ω) = 0, |ω| < ∞. Тогда H(iω) = |H(iω)|
и его импульсная характеристика фильтра равна
81
h(t) =
1
2π
+¥
ò
H(iω) eiωt dω.
-¥
Поскольку h(t) = h(t), т.е. h(t) ≠ 0 при t < 0, фильтр является физически нереализуемым (в физически реализуемом фильтре h(t) = 0
при t < 0).
Данный пример свидетельствует о том, что выбор ФЧХ формирующего фильтра должен быть сделан таким образом, чтобы фильтр
был физически реализуем. Синтез такого фильтра в общем случае
сложен. Однако, если СПМ процесса – дробно-рациональная функция вида
)
b ω2M + bM-1ω (
+  + b0
,
Sξ (ω ) = M
2(N-1)
2N
aN ω
+ aN-1ω
+  + a0
2 M-1
(8.9)
где bm , m = 0, M и an , n = 0, N – действительные числа (N ≥ M), то
процедура синтеза формирующего фильтра проста и состоит в следующем.
Из теории алгебраических уравнений известно, что неотрицательный с действительными коэффициентами полином степени 2K имеет лишь K пар комплексно сопряженных корней. Пусть
αn , α*n , n = 1, N и βm ,β*m , m = 1, M корни полиномов, стоящих в числителе и знаменателе (8.9) соответственно. Тогда (8.9) можно переписать
M
M
Õ (ω - βm ) Õ (ω - β*m )b
=1
Sξ (ω ) = mN
m=1
N
n=1
n=1
M
Õ (ω - αn ) Õ (
ω - α*n
)
aN
.
(8.10)
В силу того, что Sξ(ω) ≥ 0, aN и bM – положительные числа, и корни αn , α*n , n = 1, N и βm ,β*m , m = 1, M можно обозначить так, что αn
и βm будут иметь положительные мнимые части
(8.11)
Im {αn }³ 0, Im {βm }³ 0. Сделаем в (8.10) замену переменной ω = is в первой дроби и ω = is*
во второй (это возможно, так как ω – действительная переменная).
Тогда
2
F (s)F s*
F (s)
(8.12)
Sξ (-is) =
=
,
2
G (s)G s*
G (s)
( )
( )
82
где
F (s) =
bM M
Õ (s + fm ), aN m=1
G (s) =
N
Õ (s + gn ), (8.13)
(8.14)
n=1
fm = iβm, gn = iαn. Вследствие условия (8.14) полиномы F(s) и G(s)
имеют корни в левой полуплоскости комплексной переменной s. Поэтому, если взять в качестве частотной характеристики формирующего фильтра
H (s) =
F (s)
,
G (s)
(8.15)
мы получим устойчивый реализуемый фильтр, АЧХ которого в силу (8.12) и (8.15) соответствует условию (8.8).
Таким образом, алгоритм моделирования гауссовского случайного процесса с заданными корреляционными (спектральными)
свойствами состоит в пропускании реализации белого шума с единичной спектральной плотностью через линейный фильтр, частотная характеристика которого соответствует (8.15).
Для успешного решения задачи моделирования методом формирующего фильтра необходимо, чтобы СПМ случайного процесса на
входе описывалась дробно-рациональной функцией от переменной
ω вида (8.9). В ходе синтеза формирующего фильтра находят корни полиномов, стоящих в числителе и знаменателе (8.9), и приводят
выражения для СПМ к виду (8.10). Затем определяют корни полиномов – αn ,n = 1, N и βm ,m = 1, M, которые имеют положительные
мнимые части (см. условие (8.11)), и составляют полиномы F(s) и G(s)
в соответствии с (8.13) и (8.14). Частотная характеристика формирующего фильтра получается на основании (8.15). Необходимо заметить, что для получения отрезка случайного процесса с заданными свойствами на выходе формирующего фильтра необходимо, чтобы переходные процессы в фильтре закончились. Поэтому к сохранению выборки случайного процесса нужно перейти только после
того, как модель проработала некоторое время на «холостом ходу».
Обычно это время оценивается как tхх = (2..3)/∆fэф, где ∆fэф – эффективная ширина полосы пропускания формирующего фильтра
(см. далее).
Полученный в ходе такого синтеза фильтр является аналоговым.
Кроме того, белый шум, который подается на вход формирующе83
го фильтра, также является «аналоговым» сигналом. Поэтому необходимо создать математические модели формирующего фильтра
и белого шума. Созданию математической модели формирующего
фильтра, т.е. цифрового фильтра, посвящен следующий раздел.
Рассмотрим теперь вопрос о генерации белого шума при цифровом моделировании. Случайный процесс, который в научной литературе называется белым шумом, имеет спектральную плотность
мощности N0/2, постоянную во всей частотной области. Следовательно, его средняя мощность (дисперсия) бесконечна. Поэтому белый шум не является физически реальным процессом, а представляет собой удобную математическую абстракцию. При создании его
цифровой модели необходимо сохранить два его основных свойства:
постоянство СПМ в частотной области и статистическую независимость временных отсчетов, взятых в произвольные моменты времени. Второе свойство реализуется при моделировании весьма просто:
за реализацию дискретного белого шума берется набор независимых случайных чисел, получающихся на выходе генератора случайных нормально распределенных величин с нулевым математическим ожиданием и некоторой дисперсией σ2. Очевидно, что величина этой дисперсии (т.е. средняя мощность процесса) должна быть
каким-то образом связана со спектральной плотностью мощности
«аналогового» белого шума N0/2, дискретную модель которого мы
создаем. Для того чтобы связать эти величины, воспользуемся требованием постоянства СПМ в частотной области.
Известно, что при дискретизации непрерывного сигнала с периодом взятия отсчетов Т, спектральная функция дискретизированного сигнала становится периодической с периодом 2π/Т. Поэтому,
если СПМ нашего дискретного белого шума будет постоянна и равна N0/2 на интервале [–π/Т, π/Т], то автоматически она будет постоянна и во всей частотной области. С другой стороны, дисперсия
случайного процесса с равномерной в указанном интервале СПМ
равна (8.2)
1
σ =
2π
2
π/T
ò
-π/T
N0
N
dω = 0 .
2
2T
Отсюда при N0/2 = 1 получаем, что дисперсия независимых случайных чисел, составляющих реализацию дискретного белого шума, должна быть равна
84
σ 2 = T -1 . Следовательно, реализация дискретного белого шума размером
в N отсчетов должна быть выполнена в соответствии со следующим
равенством
1
ζ[n ] =
z[n ], n = 1, N,
T
где z[n] – случайные независимые числа с нормальным стандартным распределением (нормальное стандартное распределение – распределение Гаусса с нулевым математическим ожиданием и единичной дисперсией). Числа z[n] можно получить, используя встроенные процедуры генерации. В среде MATLAB такой процедурой
является функция randn(…).
Цифровое моделирование линейных фильтров
При цифровом моделировании линейного аналогового фильтра
его входной x(t) и выходной y(t) сигналы представляются в виде решетчатых функций x[n] и y[n], которые отличны от нуля для дискретных моментов времени nТ, n = 0,1,2,..., где Т – период дискретизации. Точное равенство входных x[n] = x(nТ) и выходных сигналов y[n] = y(nТ) аналогового фильтра и его цифровой модели не
достижимо. На практике возможно лишь приближенное равенство
y[n] ≈ y(nТ). Задачей синтеза цифровой модели является нахождение такого алгоритма вычисления y[n] по x[n], при котором это приближенное равенство выполняется как можно точнее. Рассмотрим
некоторые методы составления цифровой модели линейных аналоговых фильтров.
Метод инвариантности импульсной характеристики. При
синтезе модели этим методом обеспечивается равенство импульсных характеристик аналогового фильтра и цифрового фильтра
h[n ] = h (nT ), n = 0,1,. (8.16)
Передаточная функция цифрового фильтра при этом равна Z
-преобразованию от импульсной характеристики
H (z) =
¥
å h[n]z-n . (8.17)
n=0
Допустим, что аналоговый фильтр имеет передаточную функцию вида
H (s) = H0
(s - S1 )(s - S2 )(s - SM )
,
(s - s1 )(s - s2 )(s - sN ) 85
где Sm , m = 1, M и sn , n = 1, N – нули и полюсы передаточной характеристики. Импульсная характеристика аналогового фильтра однозначно определяется как обратное преобразование Лапласа от H(s),
и в случае, когда полюсы sn , n = 1, N простые, равна
N
h (t) = å An esnt ,
n=1
где An = H (s)(s - sn ) s=sn ; Re {sn }< 0, n = 1, N . Тогда на основании
(8.22) импульсная характеристика цифрового фильтра равна
N
h[n ] = å An esnT n . (8.18)
n=1
Подставляя (8.18) в (8.17), получим коэффициент передачи цифрового фильтра
N
An ∆t
.
H (z) = å
sn ∆t -1
z
n=1 1 - e
Учитывая, что умножение Z-преобразования цифрового сигнала на z–1 эквивалентно задержке сигнала на период дискретизации
Т, получим схему реализации фильтра в виде N параллельно включенных цифровых звеньев первого порядка.
Рассмотренный метод не может быть использован, когда коэффициент передачи аналогового фильтра H(s) не имеет полюсов (N = 0)
или когда число полюсов меньше числа нулей (N < M).
Метод билинейного преобразования. Метод используется, когда передаточная функция аналогового фильтра задана в виде
b s M + bM-1s M-1 +  + b0
H (s) = M N
aN s + aN-1s N-1 +  + a0
(8.19)
и заключается в замене
s=
2 1 - z-1
.
T 1 + z-1
(8.20)
Метод билинейного преобразования гарантирует хорошее совпадение частотных характеристик аналогового и цифрового фильтров
в диапазоне частот ω £ π 4T . При ω > π 4∆t частотные характеристики фильтров, как правило, значительно различаются.
Метод замены дифференциалов. Метод используется, когда передаточная функция аналогового фильтра задана в виде (8.19) и заключается в замене
86
1 - z-1
.
(8.21)
T
Метод замены дифференциалов применим в ситуациях, когда
метод инвариантности импульсной характеристики не используется. В то же время метод может дать значительную разницу в частотных характеристиках фильтров, если нули коэффициента передачи
аналогового фильтра H(s) имеют действительную часть, превышающую по величине 1/2T. Кроме того, метод не всегда гарантирует
близкое совпадение характеристик аналогового и цифрового фильтра, когда H(s) не имеет нулей (M = 0).
s=
Индивидуальные задания к лабораторной работе
Таблица 8.1
№
1
2
3
4
СПМ процесса
S (ω ) =
S (ω ) =
S (ω ) =
S (ω ) =
(ω
Параметры*)
k
2
+ Ω2H
)
k, ΩH
k, ΩH
k
2
(ω2 + Ω2H )
k, ΩH
k
(ω
2
3
+ Ω2H
)
k, ΩH
k
4
(ω2 + Ω2H )
5
(ω2 + Ω2B )
S (ω ) = k
2
(ω2 + Ω2H )
k, ΩH, ΩB = ΩH/√ (2)
6
ω2 + Ω2B )
(
S (ω ) = k
3
(ω2 + Ω2H )
k, ΩH, ΩB = ΩH/√ (3)
(ω2 + Ω2B )
4
(ω2 + Ω2H )
k, ΩH, ΩB = ΩH/√ (4)
7
S (ω ) = k
87
Продолжение табл. 8.1
№
СПМ процесса
8
S (ω ) = k
Параметры*)
2
ω2 + Ω2B
(
)
3
(ω2 + Ω2H )
2
k, ΩH, ΩB = ΩH/√ (2)
3
k, ΩH, ΩB = √ (3/4) ΩH
9
ω2 + Ω2B )
(
S (ω ) = k
4
(ω2 + Ω2H )
10
(ω2 + Ω2B )
S (ω ) = k
4
(ω2 + Ω2H )
11
12
13
14
15
16
17
88
S (ω ) = k
S (ω ) = k
S (ω ) = k
S (ω ) = k
S (ω ) = k
S (ω ) = k
S (ω ) =
k, ΩH, ΩB = √ (2/3) ΩH
k, ΩH
ω2
2
(ω2 + Ω2H )
k, ΩH
ω2
3
ω2 + Ω2H
(
)
k, ΩH
ω2
4
(ω2 + Ω2H )
k, ΩH
ω4
3
(ω2 + Ω2H )
k, ΩH
ω4
4
ω2 + Ω2H
(
)
k, ΩH
ω6
4
(ω2 + Ω2H )
(ω
k
4
+ Ω4H
)
k, ΩH
Окончание табл. 8.1
№
18
19
20
СПМ процесса
S (ω ) = k
S (ω ) =
Параметры*)
ω + Ω2B
ω4 + Ω4H
k, ΩH, ΩB = ΩH/√ (2)
2
(
)
k, ΩH
k
2
ω4 + Ω4H
(
S (ω ) = k
)
ω2 + Ω2B
k, ΩH, ΩB = ΩH/√ (2)
2
(ω4 + Ω4H )
2
k, ΩH, ΩB = ΩH/√ (2)
3
k, ΩH, ΩB = ΩH/√ (2)
21
(ω2 + Ω2B )
S (ω ) = k
2
(ω4 + Ω4H )
22
(ω2 + Ω2B )
S (ω ) = k
2
(ω4 + Ω4H )
23
24
25
S (ω ) = k
S (ω ) = k
S (ω ) = k
ω2
k, ΩH
ω2
k, ΩH
(ω4 + Ω4H )
2
(ω4 + Ω4H )
k, ΩH
ω4
2
ω4 + Ω4H
(
)
* Примечание. Незаданные параметры k, ΩH выбираются таким образом, чтобы
средняя мощность случайного процесса была единичной
¥
P=
1
2π
ò
Sξ (ω )dω = 1, (8.22)
-¥
а эффективная ширина спектра была равна 1 Гц
∆fýô =
1
2π Smax
¥
ò
-¥
Sξ (ω )dω =
1
= 1Ãö, Smax
(8.23)
89
где Smax – максимальное значение СПМ. Для вычисления интеграла в (8.22) и (8.23) целесообразно воспользоваться следующей формулой:
¥
ò
-¥
где
Qn (ω )
π i Mn
dω =
,
Pn (ω )Pn (-ω )
B0 Nn
Qn (ω ) = A0 ω2n-2 + A1ω2n-4 +  + An-1,
Pn (ω ) = B0 ωn + B1ωn-1 +  + Bn ,
причем все корни Pn(z) лежат в верхней полуплоскости,
A0
B0
Mn =
0
0
A1
B2
B1
0
A2
B4
B3
0
... An-1
...
0
; Nn =
...
0
... Bn
B1
B0
0
0
B3
B2
B1
0
B5
B4
B3
0
... 0
... 0
.
... 0
... Bn
В частности для n = 1,2,3,4 определители Mn и Nn равны
M1 = A0 ; N1 = B1;
M2 = A1 B0 - A0 B2 ; N2 = B1 B2 ;
M3 = B3 ( A1 B0 - A0 B2 )- A2 B0 B1; N3 = B3 (B0 B3 - B1 B2 );
M4 = A0 B4 (B2 B3 - B1 B4 )+ B0 B4 ( A2 B1 - A1 B3 )- A3 B0 (B0 B3 - B1 B2 );
(
)
N4 = B4 B0 B32 + B12 B4 - B1 B2 B3 B4
Порядок выполнения и требования к отчету
В соответствии с номером варианта необходимо сгенерировать
выборку случайного процесса с заданной спектральной плотностью
мощности двумя из рассмотренных методов. Рекомендуется взять
размер выборки (количество отсчетов генерируемого случайного
процесса) K = 512.
В качестве примера используйте файлы LAB2.M и LAB3.M (см.
далее). Рекомендуется создать копию этих файлов под своим именем в каталоге \MATLAB\WORK для модификации программы в
соответствии с заданием.
Отчет по лабораторной работе должен включать:
1. Результаты аналитического расчета параметров k, ΩH и выбора периода дискретизации Т; коэффициенты передачи аналогового
90
и цифрового формирующих фильтров с обоснованием метода синтеза и таблицу рассчитанных коэффициентов цифрового фильтра.
2. Листинг программы и графики, которые должны содержать
кривые заданной плотности мощности и ее оценки, полученной по
результатам моделирования выборки случайного процесса.
3. Выводы по работе.
Пример программы генерации случайного процесса
методом ДПФ
%
%
%
%
%
%
%
ОСНОВЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ РЭС
Листинг программы LAB8.M к лабораторной работе #8
Программа моделирования случайного комплексного
нормального процесса с корреляционной функцией вида
R(t) = exp{-[(t/tau)^2 + i*w0*t]},
где tau – постоянная времени, w0 = 2*pi*f0- несущая
частота.
% Закрываем все открытые окна вывода
close all
% Очищаем рабочую область (WORKSPACE)
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Исходные данные
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T – полное время наблюдения, с
T = 1;
% tau – постоянная времени процесса, с
tau = 0.20;
% f0 – несущая частота, Гц
f0 = 20;
% N – количество отсчетов в реализации случайного
% процесса
N = 1024;
% K – количество реализаций случайного процесса
K = 20;
% dt – шаг дискретизации по времени
dt = T/N;
% df – шаг дискретизации по частоте
df = 1/(N*dt);
% t – массив моментов взятия отсчетов случайного
% процесса
t = (0:dt:T-dt)’;
% f – массив значений частот, в которых вычисляются
% отсчеты спектра случайного процесса
f = (-0.5/dt:df:0.5/dt-df)’;
% R – массив значений корреляционной функции случайного
% процесса
91
R = exp(-(t(1:N/2 + 1)/tau).^2);
R = R.*exp(-i*2*pi*f0*t(1:N/2 + 1));
R = [R; conj(R(N/2:-1:2))];
% Sp – массив значений спектра случайного процесса
Sp = abs(ifft(R));
R = [R(N/2 + 1:N); R(1:N/2)];
% Вывод графика спектра случайного процесса
figure(1)
plot(f,fftshift(Sp))
title(‘Спектральная плотность процесса’)
xlabel(‘f, Гц’)
ylabel(‘S(f)’)
grid
% Dp – массив СКО амплитуд случайных гармоник в составе
% процесса
Dp = sqrt(Sp)/sqrt(2);
Re = zeros(N,1);Ree = Re;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Цикл генерации реализаций случайного процесса %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% M – количество случайных гармоник в составе процесса
M = N/2;
% Начало цикла генерации реализаций случайного процесса
for k = 1:K
% Ar, Ai – массивых нормальных независимых сл. чисел
Ar = randn(size(Dp(1:M)));Ai = randn(size(Dp(1:M)));
% Aleft, Aright – массивы комплексных нормальных
% независимых сл. чисел
Aleft = Ar + i*Ai;
Ar = randn(size(Dp(1:M)));Ai = randn(size(Dp(1:M)));
Aright = Ar + i*Ai;
% Aa – массив амплитуд случайных гармоник в составе
% процесса
Aa = [Aleft.*Dp(1:M); zeros(N-2*M,1);...
Aright.*Dp(N-M + 1:N)];
% S – k-я реализация случайного процесса, полученная
% с использованием алгоритма FFT
S = fft(Aa);
% Вывод графика k-й реализации случайного процесса
figure(2)
plot(t,[real(S) imag(S)]);
title([‘Реализация случайного процесса,...
k = ’ int2str(k)])
xlabel(‘t, с’)
ylabel(‘s(t)’)
grid
pause(1)
92
% Оценка корреляционной функции методом быстрой свертки
Aa = Aa.*conj(Aa);
% Ree – оценка КФ по k-й реализации случайного
% процесса
Ree = fft(Aa);
% Re – массив, в котором накапливаются оценки КФ для
% последующего усреднения
Re = Re + Ree;
end
Re = [Re(N/2 + 1:N); Re(1:N/2)];
% Вычисление среднего для оценки КФ
Re = Re./K;
% R массиву истинных значений КФ добавляется массив
% оценок
R = [R Re];
% Вывод результатов моделирования:
% – истинная КФ, синий цвет
% – оценка КФ, зеленый цвет
figure(3)
title(‘Результаты моделирования: КФ и ее оценка’)
subplot(211);plot(t-0.5*T,real(R))
xlabel(‘\tau, c’)
ylabel(‘Re\{R(\tau)\}’)
grid
legend(‘истинная КФ’,’оценка КФ’)
subplot(212);plot(t-0.5*T,imag(R))
xlabel(‘\tau, c’)
ylabel(‘Im\{R(\tau)\}’)
grid
legend(‘истинная КФ’,’оценка КФ’)
93
Пример программы генерации случайного процесса
методом формирующего фильтра
%
%
%
%
%
%
ОСНОВЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ РЭС
Листинг программы LAB9.M к лабораторной работе #9
Программа моделирования случайного действительного
нормального процесса с
корреляционной функцией (КФ) вида R(t) = exp{-|t/tau|},
где tau – постоянная времени.
% Закрываем все открытые окна вывода
close all
% Очищаем рабочую область (WORKSPACE)
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Исходные данные
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T – полное время наблюдения, с
T = 1;
% tau – постоянная времени процесса, с
tau = 0.20;
% N – количество отсчетов в реализации случайного
% процесса
N = 1024;
% K – количество реализаций случайного процесса
K = 40;
% dt – шаг дискретизации по времени
dt = T/N;
% df – шаг дискретизации по частоте
df = 1/(N*dt);
% t – массив моментов взятия отсчетов случайного
% процесса
t = (0:dt:T-dt)’;
% f – массив значений частот, в которых вычисляются
% отсчеты спектра случайного процесса
f = (-0.5/dt:df:0.5/dt-df)’;
% R – массив значений корреляционной функции случайного
% процесса
R = exp(-abs(t(1:N/2 + 1)/tau));
R = [R; conj(R(N/2:-1:2))];
% Sp – массив значений спектра случайного процесса
Sp = abs(ifft(R));
R = [R(N/2:N); R(1:N/2)];
% Re – массив, в котором накапливаются оценки КФ для
% последующего усреднения
Re = zeros(size(R));
% Вывод графика спектра случайного процесса
figure(1)
94
plot(f,fftshift(Sp))
title(‘Спектральная плотность процесса’)
xlabel(‘f, Гц’)
ylabel(‘S(f)’)
grid
% Nidle – количество отсчетов на выходе формирующего
% фильтра необходимое для
% окончания переходных процессов в нем (длительность
% холостого хода)
Nidle = round(3*tau/dt);
% W – полоса пропускания формирующего фильтра
W = 1/tau;
% a, b – коэффициенты формирующего фильтра
a = (1-0.5*W*dt)/(1 + 0.5*W*dt);
b = sqrt(0.5*W)*dt/(1 + 0.5*W*dt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Цикл генерации реализаций случайного процесса
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:K
% x – массив отсчетов белого гауссовского шума
% (входной сигнал)
x = randn(N + Nidle,1)/sqrt(dt);
% y1 – предыдущий отсчет сигнала на выходе
y1 = 0;
% Цикл холостого хода фильтра
for n = 2:Nidle
y1 = a*y1 + b*(x(n) + x(n-1));
end;
% Цикл генерации случайного процесса
y(1,1) = y1; % Задание начального условия
for n = 2:N
% Разностное уравнение фильтра
y(n,1) = a*y(n-1) + b*(x(Nidle + n) + x(Nidle + n-1));
end;
% Вывод графика k-й реализации случайного процесса
figure(2)
plot(t,real(y));
title([‘Реализация случайного процесса, k = ’ ...
int2str(k)])
xlabel(‘t, с’)
ylabel(‘s(t)’)
grid
pause(1)
% Оценка корреляционной функции методом прямой свертки
% Ree – оценка КФ по k-й реализации случайного
% процесса
95
Ree = xcorr(y,y,N/2,’unbiased’);
% Накопление оценки КФ для последующего усреднения
Re = Re + Ree;
end
% Вычисление среднего для оценки КФ
Re = Re./K;
% R массиву истинных значений КФ добавляется массив
% оценок
R = [R Re];
% Вывод результатов моделирования:
% – истинная КФ, синий цвет
% – оценка КФ, зеленый цвет
figure(3)
plot([t;T]-0.5*T,real(R))
title(‘Результаты моделирования: КФ и ее оценка’)
xlabel(‘\tau, c’)
ylabel(‘R(\tau)’)
grid
legend(‘истинная КФ’,’оценка КФ’)
Контрольные вопросы
1. Какой физический смысл имеет спектральная плотность мощности?
2. Каковы основные свойства СПМ?
3. Определите основные шаги генерации отсчетов случайного
стационарного процесса методом ДПФ.
4. Что такое алгоритм БПФ и каковы его свойства?
5. Как связаны между собой СПМ процессов на входе и выходе
линейного фильтра?
6. Каким образом моделируются дискретный белый шум? Каковы свойства этого дискретного случайного процесса?
7. Определите основные шаги генерации отсчетов случайного
стационарного процесса методом формирующего фильтра.
Литература: [7, с. 55–163; 9, c. 30–43; 14]
96
Лабораторная работа № 10
Моделирование системы автоматической
регулировки усиления
Цель работы: исследование динамических характеристик системы автоматической регулировки усиления (АРУ) путем математического моделирования.
Математическая модель системы автоматической
регулировки усиления
Система АРУ служит для согласования динамических диапазонов изменения входных и выходных сигналов усилительных каналов. Структурная схема усилителя, охваченного петлей АРУ, приведена на рис. 10.1.
Радиосигнал uвх(t) поступает на сигнальный вход усилителя, коэффициент усиления которого изменяется под воздействием сигнала Uр(t), приходящего из петли АРУ на управляющий вход. Регулирующий сигнал Uр(t) формируется на выходе петли АРУ, которая
включает: амплитудный детектор АД, схему задержки СЗ, усилитель постоянного тока УПТ и фильтр нижних частот ФНЧ.
Петля АРУ работает следующим образом. Сигнал АД UАД, величина которого пропорциональна амплитуде выходного сигнала
uвых(t), поступает на схему задержки СЗ, где производится сравнение этого сигнала с напряжением задержки Ез. Схема задержки необходима для предотвращения работы петли АРУ при малых входных сигналах. Сигнал на выходе СЗ равен
ìïUÀÄ , UÀÄ ³ Eç
UÑÇ = ïí
.
ïïî0, UÀÄ < Eç
Этот сигнал усиливается в УПТ и поступает на ФНЧ, выходным
сигналом которого является регулирующее напряжение Uр(t). МожV»ÎU
V»ÔÎU
¬ÊÁÄÁ˾ÄÕ
VQU
­¦°
™
¬¨«
ª &À
Рис. 10.1. Структурная схема усилителя,
охваченного петлей АРУ
97
но считать, что быстродействие системы АРУ определяется характеристиками ФНЧ.
Исходными данными для расчета системы АРУ являются:
– минимальное Uвх min и максимальное Uвх max значения амплитуды входного сигнала;
– минимальное Uвых min и максимальное Uвых max значения амплитуды выходного сигнала.
На основании данных величин можно вычислить динамические
диапазоны входных и выходных сигналов
Dâõ = 20 lg
Uâûõ max
Uâõ min
, Dâûõ = 20 lg
Uâûõ max
Uâûõ min
.
Учитывая, что максимальный и минимальный коэффициенты
усиления равны
Kmax = 20 lg
Uâûõ min
Uâõ min
, Kmin = 20 lg
Uâûõ max
Uâõ max
,
несложно показать, что
Dâõ - Dâûõ = Kmax - Kmin . Отсюда, зная Kmax, Dвх и Dвых, можно определить Kmin.
Для анализа и моделирования системы АРУ необходимо знать
две характеристики: регулировочную и амплитудную. Регулировочной характеристикой (РХ) называется зависимость коэффициента усиления K от регулирующего напряжения Uр. Примерный
вид РХ представлен на рис. 10.2.
Как правило, зависимость K = K(Up) имеет близкий к линейному
характер (рис. 10.2, а). Поэтому простейшей аппроксимацией РХ,
пригодной для моделирования, является линейная
¹
º
,
, NBY
6Q
6QNBY
, NJO
6Q
G
6QNBY
6»ÔÎ
6»ÔÎ NJO
6»ÔÎ NBY
Рис. 10.2. Регулировочная и амплитудная характеристики АРУ
98
K = Kmax - µ U p
где µ = tgγ – крутизна РХ. Возможны и другие аппроксимации: экспоненциальная, полиномиальная и т.д. По РХ определяется максимальное регулирующее напряжение Up max.
Зависимость Up = Up(Uвых) называется амплитудной характеристикой (АХ) (рис. 10.2, б). Эта характеристика также имеет близкий
к линейному характер. Статический коэффициент усиления петли
АРУ определяется по АХ как
KÀÐÓ =
Up max
Uâûõ max - Uâûõ min
=
(
Up max
.
)
Eç 100,05Dâûõ -1
С другой стороны
K
=K
K
,
ÀÐÓ
ÀÄ ÓÏÒ где KАД = 0.8…0.9 – коэффициент усиления АД; KУПТ – коэффициент усиления УПТ. Отсюда можно найти потребный коэффициент
усиления УПТ
KÓÏÒ =
Up max
(
.
)
Eç KÀÄ 100,05Dâûõ -1
Отметим, что из последнего уравнения следует вывод о невозможности построения идеальной АРУ с Dвых = 0 дБ, так как для
этого необходим усилитель постоянного тока с бесконечно большим
коэффициентом усиления.
Выбор ФНЧ в петле АРУ влияет не только на быстродействие системы, но и на ее устойчивость. В простейшем случае ФНЧ представляет собой интегрирующую RC-цепочку. Однако возможно использование и более сложных фильтров.
Анализ динамического режима работы АРУ является сложной
задачей, поскольку данная система представляет собой замкнутую
нелинейную систему регулирования. Полностью провести такой
анализ возможно лишь в некоторых простейших ситуациях. Решим, например, такую задачу: требуется определить реакцию системы АРУ на скачок амплитуды входного сигнала при условии,
что регулировочная характеристика линейна, а ФНЧ представляет
собой фильтр первого порядка с постоянной времени τ (т.е. ФНЧ –
интегрирующая RC-цепочка).
99
Найдем дифференциальное уравнение (ДУ), описывающее работу системы АРУ в рассматриваемом случае. Начнем с ДУ, связывающего сигналы на входе и выходе ФНЧ:
dUp (t)
+ Up (t) = KÀÐÓ éëUâûõ (t)- Eç ùû ,
dt
где учтено, что входным сигналом фильтра является разность напряжений [Uвых(t) – Ез], прошедшая через усилительное звено с коэффициентом усиления KАРУ. С другой стороны, амплитуды выходного и входного сигналов связаны между собой следующим уравнением:
τ
Uâûõ (t) = K (t)Uâõ (t) = éêë Kmax - µUp (t)ùúû Uâõ (t).
Данное уравнение можно решить относительно регулирующего
напряжения Uр(t)
U
(t)ù 1
1é
Up (t) = êê Kmax - âûõ úú = éë Kmax - K (t)ùû ,
µ ëê
Uâõ (t) ûú µ
Подставляя данное выражение в ДУ, получим
τ
dK (t)
dt
+ éë1 + µKÀÐÓ Uâõ (t)ùû K (t) = Kmax + µKÀÐÓ Eç .
Учитывая, что входной сигнал в нашем случае равен скачку напряжения
ì
ïUmax , t ³ 0
Uâõ (t) = ï
,
í
ï
0, t < 0
ï
î
окончательно получим следующее ДУ для коэффициента усиления
τ
dK (t)
dt
+ [1 + µKÀÐÓ Umax ]K (t) = Kmax + µKÀÐÓ Eç .
Данное уравнение является линейным ДУ с постоянными коэффициентами, и его решение может быть легко вычислено. Учитывая, что в начальный момент времени K(0) = Kmax, в результате вычислений имеем
Kmax + µKÀÐÓ Eç
1 - e-t/T ,
1 + µKÀÐÓ Umax
τ
T=
.
1 + µKÀÐÓ Umax
K (t) = Kmax e-t/T +
100
(
)
6»Î<L>
6»ÔÎ <L>
,<L>
s
&À
[
6 <L>
Q
¨­
)[
,™©¬
Рис. 10.3. Функциональная схема математической
модели усилителя, охваченного петлей АРУ
Параметр T представляет собой по сути постоянную времени
АРУ, которая, как следует из предыдущего уравнения, меньше
постоянной времени фильтра, причем это уменьшение зависит от
уровня входного сигнала Umax: чем больше сигнал, тем меньше постоянная времени Т, и, следовательно, больше быстродействие системы АРУ.
В ситуациях, когда РХ или АХ являются существенно нелинейными функциями или ФНЧ имеет высокий порядок, аналитические расчеты динамических параметров АРУ представляют большую сложность. Поэтому единственным средством для их оценки
является математическое моделирование. Построение математической модели системы АРУ не вызывает затруднений. Функциональная схема математической модели усилителя, охваченного петлей
АРУ, приведена на рис. 10.3.
Входным сигналом АРУ является амплитуда входного сигнала
Uвх[k], выходным сигналом – амплитуда выходного сигнала усилителя Uвых[k], связь между которыми имеет вид
Uâûõ [k] = K [k]Uâõ [k],
где K[k] – значение коэффициента усиления. Коэффициент усиления
получается, как результат работы петли АРУ, включающей: вычитатель, на второй вход которого подается напряжение задержки Ез;
схему задержки; усилитель с коэффициентом усиления KАРУ; ЦФ с
коэффициентом передачи H(z); нелинейное безынерционное звено с
характеристикой, соответствующей РХ АРУ; элемент временной задержки на один период дискретизации. Схема задержки представляет собой нелинейное безынерционное звено с характеристикой
ìïuâõ ÑÇ , uâõ ÑÇ ³ 0
uâûõ ÑÇ = g uâõ ÑÇ = ïí
,
ïï0, uâõ ÑÇ < 0
î
(
)
101
где uвх СЗ и uвых СЗ – входной и выходной сигналы СЗ. Звено РХ также является нелинейным и безынерционным. Для его реализации
можно использовать либо реальную РХ моделируемой системы АРУ,
либо ее аппроксимацию (линейную, экспоненциальную, полиномиальную и др.). Звено временной задержки на один такт необходимо
для упрощения вычислений при моделировании, поскольку в данном случае отпадает необходимость решать нелинейное уравнение
относительно текущего значения коэффициента усиления [9, с. 58].
Цифровой фильтр является моделью ФНЧ петли АРУ, и ее построение не вызывает затруднений.
Индивидуальные задания к лабораторной работе
Таблица 10.1
№
1
2
3
4
Задание
Провести сравнительный анализ результатов моделирования
линейной системы АРУ первого порядка с результатами теоретического расчета
Исследовать влияние нелинейности регулировочной характеристики на переходные процессы в системе АРУ
Исследовать переходные процессы в АРУ второго порядка
Исследовать переходные процессы в АРУ третьего порядка
Порядок выполнения и требования к отчету
В соответствии с номером задания необходимо:
– составить план проведения исследований и согласовать его с
преподавателем;
– провести необходимые расчеты параметров АРУ при следующих исходных данных: динамический диапазон входных сигналов
60 дБ, динамический диапазон выходных сигналов 20 дБ, максимальный коэффициент усиления 20 дБ, максимальное регулирующее напряжение 2 В, напряжение задержки 0.01 В, постоянная времени ФНЧ 0.01 с.
В качестве примера используйте файл LAB10.M (см. далее). Рекомендуется создать копию этого файла под своим именем в каталоге \
MATLAB\WORK для модификации программы в соответствии с заданием.
Отчет по лабораторной работе должен включать:
1. Результаты аналитического расчета параметров АРУ и выбора периода дискретизации Т, коэффициенты передачи цифровых
фильтров.
102
2. Листинг программы и графики, содержащие результаты моделирования.
3. Теоретические расчеты зависимостей, исследованных в ходе
выполнения работы.
4. Выводы по работе.
Пример программы моделирования системы АРУ
%
%
%
%
%
ОСНОВЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ РЭС
Листинг программы к лабораторной работе #10
Программа моделирования системы АРУ:
ФНЧ – 1-го порядка,
РХ и АХ – линейные функции.
% Закрываем все открытые окна вывода
close all
% Очищаем рабочую область (WORKSPACE)
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Исходные данные
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Период дискретизации, с
Ts = 1e-3;
% Динамический диапазон входных сигналов, дБ
DindB = 60;
% Динамический диапазон выходных сигналов, дБ
DoutdB = 20;
% Максимальный коэффициент усиления, дБ
KmaxdB = 20;
% Максимальное регулирующее напряжение, В
Upmax = 2;
% Напряжение задержки, В
Ed = 0.01;
% Постоянная времени ФНЧ, с
tau = 0.01;
% Величина скачка напряжения на входе, В
Umax = 0.1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Расчет АРУ в статическом режиме работы
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Минимальный коэффициент усиления, дБ
KmindB = KmaxdB-(DindB-DoutdB);
Kmax = 10^(0.1*KmaxdB);
Kmin = 10^(0.1*KmindB);
mu = (Kmax-Kmin)/Upmax; % Крутизна РХ
% коэффициент усиления петли АРУ
Kagc = Upmax/(10^(0.05*DoutdB)-1);
% Коэффициенты ФНЧ
103
a = (1-2*tau/Ts)/(1 + 2*tau/Ts);
b = 1/(1 + 2*tau/Ts);
t = (0:1e2)’*Ts; % Текущее время
N = length(t); % Размер выборки
Uin = Umax*ones(N,1); % Входной сигнал
K = zeros(N,1); % Коэффициент усиления
K(1) = Kmax;
Uout = zeros(N,1); % Выходной сигнал
Uout(1) = Kmax*Umax;
Up = zeros(N,1); % Напряжение регулировки
uin1 = 0; % Предыдущее значение сигнала на входе ФНЧ
uout1 = 0; % Предыдущее значение сигнала на выходе ФНЧ
for n = 2:N
% Сигнал на выходе схемы задержки
uin0 = Kagc*max([0,Uout(n-1)-Ed]);
% Текущее значение сигнала на выходе ФНЧ
uout0 = -a*uout1 + b*(uin0 + uin1);
uin1 = uin0;
uout1 = uout0;
% Ограничиваем регулирующий сигнал
if uout0<0
Up(n) = 0;
elseif uout0 > = Upmax
Up(n) = Upmax;
else
Up(n) = uout0;
end
K(n) = Kmax-mu*Up(n); % Текущий КУ
% Текущее значение выходного сигнала
Uout(n) = K(n)*Uin(n);
end
% Теоретический расчет коэффициента усиления
T = tau/(1 + mu*Kagc*Umax); % Постоянная времени АРУ
% Установившееся значение КУ
C = (Kmax + mu*Kagc*Ed)/(1 + mu*Kagc*Umax);
Kth = Kmax*exp(-t/T) + C*(1-exp(-t/T));
figure(1)
plot(t,K,t,Kth)
grid
xlabel(‘t’)
ylabel(‘K(t)’)
title(‘Коэффициент усиления’)
legend(‘модель’,’теория’)
figure(2)
plot(t,Up)
grid
104
xlabel(‘t’)
ylabel(‚U_p(t)‘)
title(‚Регулирующее напряжение’)
figure(3)
plot(t,Uout)
grid
xlabel(‘t’)
ylabel(‘U_{вых}(t)’)
title(‘Выходной сигнал’)
Контрольные вопросы
1. Каково назначение системы АРУ?
2. Нарисуйте функциональную схему системы АРУ и определите назначение основных элементов.
3. Что такое регулировочная и амплитудная характеристики
АРУ и какой вид они имеют?
4. Нарисуйте функциональную схему математической модели
системы АРУ и определите назначение основных элементов.
5. Как рассчитываются параметры математической модели системы АРУ в статическом режиме работы?
Литература: [9, с. 82–86; 15; 16, с. 158–178]
105
Лабораторная работа № 11
Моделирование системы фазовой
автоматической подстройки частоты
Цель работы: исследование характеристик системы фазовой автоматической подстройки частоты (ФАПЧ) путем математического
моделирования.
Математическая модель системы фазовой
автоматической подстройки частоты
Система ФАПЧ служит для обеспечения равенства частот сигнала, поступающего на вход системы, и сигнала, генерируемого
управляемым генератором. Подобная задача встречается в различных ситуациях: подстройка частоты гетеродина преобразователя
частот, выделение (измерение) частоты принимаемых сигналов,
синхронизация телекоммуникационных систем и др. Структурная
схема ФАПЧ приведена на рис. 11.1.
Система ФАПЧ представляет собой замкнутую нелинейную систему автоматического слежения за фазой входного сигнала. Основными элементами ФАПЧ являются фазовый детектор (ФД), фильтр
низкой частоты (ФНЧ), управляющий элемент (УЭ) и управляемый
генератор (УГ). На вход ФД поступает сигнал uc(t) = Uccos(ω0t + j0),
на второй вход ФД подается сигнал УГ uг(t) = Uгsin(ωt + j). На выходе ФНЧ выделяется сигнал равный uу(t) = 0.5sin(ψc – ψг), где
ψc = ω0t + j0 – фаза входного сигнала, ψг = ωt + j – фаза сигнала
УГ. Данный сигнал поступает на вход УЭ, воздействие которого непосредственно на колебательную систему УГ приводит к изменению
частоты его сигнала. С течением времени между фазами входного
сигнала и сигнала УГ устанавливается постоянная разность, и, следовательно, обеспечивается равенство частот ω0 и ω.
Для синтеза математической модели ФАПЧ линеаризируем ее и
рассчитаем основные параметры. Будем считать, что разность фаз
(ψc – ψг) на столько мала, что можно без особой погрешности полоVD U
­
V¼ U
¬œ
­¦°
VÌU
¦
V¼ U
Рис. 11.1. Структурная схема ФАПЧ
106
YDU
)T
Y¼U
MT
Рис. 11.2. Линеаризированная функциональная схема ФАПЧ
жить sin(ψc – ψг) ≈ (ψc – ψг), а частота УГ линейно зависит от управляющего напряжения. Тогда функциональная схема линеаризованной системы ФАПЧ будет иметь вид, представленный на рис. 11.2.
В этой схеме H(s) – коэффициент передачи ФНЧ, µ – коэффициент преобразования управляющего напряжения в фазу УГ, 1/s – коэффициент передачи интегратора, наличие которого в схеме соответствует интегральной зависимости фазы сигнала УГ от значений
его мгновенной частоты. Будем считать, что ФНЧ является фильтром первого порядка с постоянной времени τ. Тогда коэффициент
передачи разомкнутого контура ФАПЧ будет равен
K0 (s) =
µ
1
,
2 s(τs + 1)
(11.1)
а коэффициент передачи замкнутого контура
K (s) =
K0 (s)
=
Ω20
µ
,
= 2
2s(τs + 1)+ µ s + 2ξΩ0 s + Ω20
(11.2)
1 + K0 (s)
где введены следующие обозначения: Ω0 = µ/2τ и 2ξΩ0 = 1/τ. Введение
данных обозначений позволило привести запись коэффициента передачи ФАПЧ (11.1) к виду, принятому для колебательного звена.
При этом введенные параметры имеют следующий смысл: Ω0 – граничная частота полосы пропускания звена; ξ – коэффициент переколебательности.
Известно, что переходная характеристика колебательного звена
равна
é
ù
ξ
ê
ú
g (t) = 1 - e-ξΩ0t ê cos Ω0 1 - ξ2 t +
sin Ω0 1 - ξ2 t ú .
ê
ú
1 - ξ2
êë
úû (
)
(
)
Примерный вид переходной характеристики приведен на рис. 11.3.
Величина первого выброса g(t) над единичным уровнем называется перерегулированием. Можно показать, что
é
ù
πξ ú
ê
∆ = exp ê(11.3)
ú. ê 1 - ξ2 ú
êë
úû
107
HU
$
U
Рис. 11.3. Переходная характеристика
колебательного звена
При выборе динамических характеристик замкнутой системы
слежения задают граничную частоту полосы пропускания Ω0 и перерегулирование ∆. При этом коэффициент переколебательности
может быть вычислен на основании (11.3) как
ξ=
1
2
1 + (π ln ∆ )
.
(11.4)
Обычно ∆ = 0.3…0.4. При этом ξ = 0.28…0.36. Зная Ω0 и ξ, могут
быть вычислены параметры:
τ=
1
1
и µ= .
ξ
2ξΩ0
(11.5)
На этом выбор параметров математической модели ФАПЧ можно считать завершенным и можно составить ее функциональную
схему, которая представлена на рис. 11.4.
Данная функциональная схема содержит:
– нелинейный дискриминатор, который вычисляет sin(ψc – ψг),
– линейный цифровой фильтр, коэффициент передачи которого
K0(z) имеет прототипом аналоговый фильтр с коэффициентом передачи K0(s) (см. уравнение (11.3)),
– вычитатель и элемент задержки, который включен в схему для
упрощения вычислений при моделировании (см. [9, с. 58]).
YD<O>
TJOY D <O>sY ¼ <Os>
,[
[s
Рис. 11.4. Функциональная схема
математической модели ФАПЧ
108
Y¼<O>
Для определения K0(z) по соответствующему коэффициенту передачи аналогового прототипа K0(s) может быть использован метод
билинейного преобразования [6, с. 46; 9, c. 48].
Для повышения быстродействия системы ФАПЧ ФНЧ фильтр
заменяют на пропорционально-интегрирующий с коэффициентом
передачи
τ s +1
H (s) = 1
,
τs + 1 где τ1 – постоянная времени, которая выбирается для обеспечения
нужного времени захвата.
Индивидуальные задания к лабораторной работе
Таблица 11.1
№
Задание
1
Исследовать переходный процесс в ФАПЧ при захвате частоты сигнала, определить полосу и время захвата.
Исследовать переходный процесс в ФАПЧ при удержании частоты
сигнала, определить полосу удержания.
Оценить искажения формы выходного сигнала при воздействии на
ФАПЧ фазомодулированного сигнала.
Исследовать переходный процесс в ФАПЧ с пропорциональноинтегрирующим фильтром при захвате частоты сигнала, определить
полосу и время захвата.
Исследовать переходный процесс в ФАПЧ с пропорциональноинтегрирующим фильтром при удержании частоты сигнала, определить полосу удержания.
Определить статистические характеристики (плотность распределения, математическое ожидание и дисперсию) разности фаз входного
сигнала и сигнала УГ (ψc – ψг) при наличии на входе белого шума.
2
3
4
5
6
Порядок выполнения и требования к отчету
В соответствии с номером задания необходимо:
– составить план проведения исследований и согласовать его с
преподавателем;
– провести необходимые аналитические расчеты параметров
ФАПЧ при следующих исходных данных: частота входного сигнала f0 = 1 кГц, пределы изменения частоты ±50 Гц, время захвата не
более 20 мс;
– написать собственную программу моделирования в соответствии с заданием и разработанным планом исследований.
109
В качестве примера используйте файл LAB11.M (см. далее). Рекомендуется создать копию этого файла под своим именем в каталоге \
MATLAB\WORK для модификации программы в соответствии с заданием.
Отчет по лабораторной работе должен включать:
1. Результаты аналитического расчета параметров ФАПЧ и выбора периода дискретизации Т, коэффициенты передачи цифровых
фильтров.
2. Листинг программы и графики, содержащие результаты моделирования.
3. Теоретические расчеты зависимостей, исследованных в ходе
выполнения работы.
4. Выводы по работе.
Пример программы моделирования системы ФАПЧ
% ОСНОВЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ РЭС
% Листинг программы к лабораторной работе #11
% Программа моделирования системы ФАПЧ:
% ФНЧ – 1-го порядка,
% РХ и – линейная функция.
% Закрываем все открытые окна вывода
close all
% Очищаем рабочую область (WORKSPACE)
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Исходные данные
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Период дискретизации, с
Ts = 1/8e3;
% Частота сигнала, Гц
f0 = 1000;
% Начальная расстройка, Гц
df = 0;
% Полоса ФНЧ, Гц
DF = 100;
% Коэффициент перерегулирования
Delta = 0.3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Расчет параметров ФАПЧ
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xi = 1/sqrt(1 + (pi/log(Delta))^2);
W0 = 2*pi*DF;
tau = 0.5/(xi*W0);
mu = W0/xi;
% Коэффициенты ФНЧ
[numd,dend] = bilinear([0 0 1],[tau 1 0],1/Ts);
110
t = (0:2e2)’*Ts; % Текущее время
N = length(t); % Размер выборки
f = zeros(N + 2,1);
f(1) = (f0 + df);
phi_in = 2*pi*f0*t;
Uin = cos(phi_in); % Входной сигнал
phi_out = zeros(N + 2,1);
phi_out(1:3) = 2*pi*(f(1)*t(1:3));
% Предыдущее значение сигнала на входе ФНЧ
uin = sin(phi_in(4:-1:2)-phi_out(3:-1:1));
% Предыдущее значение сигнала на выходе ФНЧ
uout = zeros(2,1);
for n = 3:N-1
y = numd*uin-dend(2:3)*uout;
uout = [y;uout(1)];
phi_out(n) = 2*pi*f(1)*t(n) + mu*y;
uin = sin(phi_in(n + 1:-1:n-1)-phi_out(n:-1:n-2));
end
Uout = cos(phi_out); % Выходной сигнал
figure(1)
plot(t,phi_in,t(1:N-1),phi_out(1:N-1))
grid
xlabel(‘t’)
ylabel(‘\psi_с(t), \psi_г(t)’)
legend(‘\psi_с(t)’,’\psi_г(t)’)
figure(2)
plot(t(1:N-1),phi_in(1:N-1)-phi_out(1:N-1))
grid
xlabel(‘t’)
ylabel(‘\Delta\psi(t)’)
figure(3)
subplot(211)
plot(t,Uin)
xlim([0 0.015])
xlabel(‘t’)
ylabel(‘u_c(t)’)
subplot(212)
plot(t(1:N-1),Uout(1:N-1))
xlim([0 0.015])
xlabel(‘t’)
ylabel(‘u_г(t)’)
111
Контрольные вопросы
1. Каково назначение системы ФАПЧ?
2. Нарисуйте функциональную схему системы ФАПЧ и определите назначение основных элементов.
3. Нарисуйте линеаризированную функциональную схему системы ФАПЧ и определите коэффициенты передачи ее звеньев.
4. Нарисуйте функциональную схема математической модели
системы ФАПЧ и определите назначение основных элементов.
5. Как рассчитываются параметры математической модели контура автоматического регулирования системы ФАПЧ?
Литература: [16, с. 195–209; 17, 290–295]
Лабораторная работа № 12
Моделирование следящего моноимпульсного
амплитудного суммарно-разностного пеленгатора
Цель работы: исследование характеристик моноимпульсной системы автоматического сопровождения по углу путем математического моделирования.
Математическая модель следящего моноимпульсного
амплитудного суммарно-разностного пеленгатора
Следящий моноимпульсный амплитудный суммарно-разностный
пеленгатор (МАСРП) предназначен для слежения за радиолокационными целями по угловым координатам (азимуту и углу места). В
целях упрощения материала рассмотрим структурную схему следящего пеленгатора, осуществляющего сопровождение цели по азимуту (рис. 12.1).
Входным сигналом для системы слежения является азимутальная координата цели αц(t). Задачей системы слежения является поддержание примерного равенства
α à (t) » α ö (t), (12.1)
где αа(t) – угловая координата равносигнального направления пеленгатора. Для выполнения равенства (12.1) в пеленгаторе используются два звена: дискриминатор и сглаживающее звено. Дискриминатор служит для оценки углового рассогласования ε(t) = αц(t) – αа(t),
которое формируется вне пеленгатора и возникает естественным образом (а не как результат вычислений) вследствие неравенства αц(t)
и αа(t). Для оценки ε(t) используются радиосигналы, рассеянные
целью и принятые приемником пеленгатора. Помимо сигналов цели на входе приемника существуют помеховые сигналы (например,
внутренние шумы приемника), которые делают невозможным точное равенство углового рассогласования ε(t) и его оценки εˆ (t), коªÁ¼Æ¹ÄÔÁÈÇžÎÁ
EU
AÏU
ÁÊÃÉÁÅÁƹËÇÉ
?
EU
)T
A¸U
A¸U
Рис.
12.1. Структурная схема следящего пеленгатора
113
3
©ª¦
GA
3s$
ÈɾǺɹ
ÀÇ»¹Ë¾ÄÕ
GA
$
ªÅ
¬¨°3
œ¾Ë
™©¬
ªÅ
¬¨°$
­
Рис. 12.2. Функциональная схема дискриминатора
торая формируется как реальный электрический сигнал на выходе дискриминатора. Оценка εˆ (t) подается на сглаживающее звено,
задачей которого является фильтрация εˆ (t) и придание замкнутой
системе слежения нужного порядка астатизма. Фактическим выходным процессом системы сопровождения является азимутальное
положение равносигнального направления αа(t).
Описанная структурная схема соответствует, в принципе, любому следящему пеленгатору. Отличительной особенностью МАСРП
является то, что дискриминатор в его составе имеет функциональную схему, представленную на рис. 12.2.
Пусть сигнал от точечной цели, имеющей угловую координату
αц(t), равен eц(t)eiωt. Сигналы в парциальных каналах 1 и 2, которые имеют антенны с совмещенными фазовыми центрами, равны
e1,2 (t) = f1,2 (α ö )eö (t)eiωt , где f1,2(α) = f(α ± αсм) – диаграммы направленности (ДН) парциальных каналов, которые равны смещенной на αсм вправо и влево от
равносигнального направления ДН f(α) (рис. 12.3).
После суммарно-разностного преобразования образуются суммарный и разностный сигналы:
eΣ,∆ (t) =
1 é
iωt
ù
êë f1 (α ö )± f2 (α ö )úû eö (t)e . 2
(12.2)
Коэффициент 1 2 в (12.2) учитывает баланс мощностей в
суммарно-разностном преобразователе. Последнее уравнение удобно переписать в виде
eΣ,∆ (t) = fΣ,∆ (α ö )eö (t)eiωt ,
1 é
ù
где fΣ,∆ (α ) =
f1 (α )± f2 (α )û – диаграммы направленности суммар2ë
ного и разностного каналов (рис. 12.3).
114
GA A
G A A
G3A
A
G$A
A
Рис. 12.3. Диаграммы направленности парциальных каналов,
суммарного и разностного каналов
После преобразования на промежуточную частоту ωпр в смесителях См, на вторые входы которых подается сигнал гетеродина
Гет, сигналы поступают на усилители промежуточной частоты
УПЧ-Σ и УПЧ-∆. Особенностью этих усилителей является то, что
на их управляющие входы поступает один и тот же сигнал автоматической регулировки усиления (АРУ) с выхода УПЧ-Σ. Таким
образом происходит нормировка сигнала разностного канала по
амплитуде сигнала суммарного канала. Если считать, что АРУ
является быстродействующей, т.е. ее постоянная времени значительно меньше времени корреляции амплитудных флюктуаций
сигнала цели, то сигналы на выходе УПЧ-Σ и УПЧ-∆ можно записать в виде
eΣ,∆ (t) =
fΣ,∆ (α ö )eö (t)e
iωïðt
fΣ (α ö )eö (t)
.
Если систему АРУ нельзя считать быстродействующей, то необходимо построить ее математическую модель (см. лабораторную работу № 10).
115
Сигнал ошибки пеленгатора получается на выходе фазового детектора ФД
{ }
εˆ (t) = EΣ (t)E∆ (t)cos éëϕ Σ (t)- ϕ ∆ (t)ùû = Re eΣ* e∆ =
f∆ (α ö )
,
fΣ (α ö )
где EΣ, E∆, jΣ, j∆ – соответственно амплитуды и фазы сигналов суммарного и разностного каналов. Отношение
Π (α ) =
f∆ (α )
fΣ (α )
называется пеленгационной характеристикой (ПХ). Особенностью
ПХ является то, что это нечетная функция с линейным участком в
окрестности α = 0 (рис. 12.4). Тангенс
¨ A угла наклона ПХ в окрестности точки
α = 0 называется крутизной
G
µ = tgγ = Π ¢(0).
A
Крутизна ПХ является исключительно важной характеристикой, поскольку определяет точность измереРис. 12.4. Пеленгационная
ния угловых координат моноимпульсхарактеристика
ным пеленгатором: чем больше крутизна, тем выше точность. Если отклонение цели от равносигнального направления мало, то можно считать, что
εˆ (t) » µ éë α ö (t)- α a (t)ùû = µε (t). (12.3)
Таким образом, при пеленгации точечной цели и отсутствии шумов и помех в каналах пеленгатора сигнал ошибки с точностью до
постоянного параметра µ равен отклонению цели от равносигнального направления (РСН).
Рассмотрим теперь сглаживающее звено. Как правило, сглаживающее звено (СЗ)– это линейный фильтр, в состав которого входит
один или несколько интеграторов. Количество интеграторов определяет порядок астатизма замкнутой системы слежения. Рассмотрим простейший случай, когда замкнутая система слежения имеет память по положению и, следовательно, ее порядок астатизма равен единице. В этом случае СЗ представляет собой последовательное
соединение интегрирующего и апериодического звеньев и его коэффициент передачи равен
116
EU
A ÏU
M
˜
EU
)T
A ¹U
Рис. 12.5. Линеаризованная
функциональная схема пеленгатора
H (s) =
k
s(1 + τ s)
(12.4)
где k – статический коэффициент передачи; τ – постоянная времени
апериодического звена. Параметры k и τ выбираются таким образом, чтобы система слежения в целом обладала необходимыми динамическими характеристиками. Для определения динамических
характеристик МАСРП составим линеаризованную функциональную схему пеленгатора (рис. 12.5).
Здесь нелинейное звено – дискриминатор, в соответствии с (12.3)
заменено линейным звеном – усилителем с коэффициентом передачи µ. Учитывая, что разомкнутая система имеет коэффициент передачи
µk
K0 (s) = µH (s) =
,
s(1 + τ s)
сквозной коэффициент передачи будет равен
K (s) =
K0 (s)
Ω20
=
,
1 + K0 (s) s2 + 2ξΩ0 s + Ω20
где Ω0 = µk τ , ξ = 0,5 µ kτ . Таким образом, замкнутая система
слежения будет эквивалентна колебательному звену с граничной
частотой полосы пропускания Ω0 и коэффициентом переколебательности ξ. Переходная характеристика этого звена имеет вид
é
ù
ξ
ê
ú
g (t) = 1 - e-ξΩ0t ê cos Ω0 1 - ξ2 t +
sin Ω0 1 - ξ2 t ú .
2
ê
ú
1- ξ
ëê
ûú (
)
(
)
Примерный график g(t) изображен на рис. 12.6.
Величина первого выброса g(t) над единичным уровнем называется перерегулированием. Несложно показать, что
ü
ïìï
πξ ïïï
∆ = exp ïí(12.5)
ý. ïï
2ï
ï
1
ξ
îï
þï
117
HU
$
U
Рис. 12.6. Переходная характеристика
При выборе динамических характеристик замкнутой системы
слежения задают граничную частоту полосы пропускания Ω0 и перерегулирование ∆. При этом коэффициент переколебательности
может быть вычислен на основании (12.5) как
1
ξ=
.
2
1 + (π ln ∆ )
Обычно ∆ = 0.3…0.4, при этом ξ = 0.28…0.36. Зная Ω0 и ξ, параметры сглаживающего звена k и τ могут быть вычислены как
Ω
1
k= 0 , τ=
.
2ξµ
2ξΩ0
На этом этап выбора параметров математической модели можно
считать завершенным и следует составить функциональную схему
математической модели МАСРП (рис. 12.7).
Входным сигналом системы является дискретный процесс αц[k].
Сигнал углового рассогласования ε[k] поступает на два нелинейных
безынерционных звена, функции передачи которых соответствуют
диаграммам направленности суммарного fс и разностного fp каналов. Выходные сигналы этих звеньев поступают на два перемножителя, на вторые входы которых подается дискретный сигнал цели eц[k]. На выходах перемножителей получаются сигналы цели в
суммарном и разностном каналах eс[k], eр[k]. Эти сигналы подаются
118
AÏ<L>
AB<L>
E<L> G 
Q
FÏ<L>
FÉ<L>
3F\^
OÉ<L>
]] Gʐ
FÊ<L>
OÊ<L>
[
)[
?
E<L>
Рис. 12.7. Функциональная схема математической модели МАСРП
на сумматоры, где к сигналам цели прибавляются сигналы nс[k] и
nр[k], являющиеся дискретными моделями помех и шумов. Далее
2
вычисляются значения процессов Re ec* [k]e p [k] и ec [k] , которые
поступают на делитель. Таким образом организуется модель АРУ
и ФД. Дискретный сигнал εˆ [k] поступает на цифровой линейный
фильтр, прототипом которого является сглаживающее звено с коэффициентом передачи (12.4). Выходной сигнал цифрового фильтра
задерживается на один период дискретизации и подается на второй
вход разностного звена. Необходимость постановки дополнительного звена задержки z–1 вызвана желанием упростить вычислений
при моделировании [9, с. 58].
За пределами схемы на рис. 12.7 остались модели формирования
сигнала цели eц[k] и шумов nс[k], nр[k]. Эти сигналы моделируются
как случайные дискретные процессы с характеристиками, соответствующими статистике сигналов цели и помех.
{
}
Индивидуальные задания к лабораторной работе
Таблица 12.1
№
Задание
Исходные данные
1
Определить статистические характеристики ошибки сопровождения
цели при различных значениях
граничной ширины контура АС
2
Определить статистические характеристики ошибки сопровождения
цели при различных значениях
отношения сигнал/шум
Пеленгатор: МАСРП
ОСШ: 20 дБ
Порядок астатизма: 1
Цель: точечная, флюкт.,
∆F = 20 Гц
Пеленгатор: МАСРП
f0 = 2 Гц
Порядок астатизма: 1
Цель: точечная, флюкт.,
∆F = 20 Гц
119
Окончание табл. 12.1
№
Задание
Исходные данные
3
Определить статистические характеристики ошибки сопровождения
цели при различных значениях
угловой скорости движения цели
4
Определить статистические характеристики ошибки сопровождения
цели при различных значениях
граничной ширины контура АС
5
Определить статистические характеристики ошибки сопровождения
цели при различных значениях
отношения сигнал/шум
6
Определить статистические характеристики ошибки сопровождения
цели при различных значениях
угловой скорости движения цели
Пеленгатор: МАСРП
f0 = 2 Гц
Порядок астатизма: 1
Цель: точечная, флюкт.,
∆F = 20 Гц
Пеленгатор: МАСРП
ОСШ: 20 дБ
Порядок астатизма: 2
Цель: точечная, флюкт.,
∆F = 20 Гц
Пеленгатор: МАСРП
f0 = 2 Гц
Порядок астатизма: 2
Цель: точечная, флюкт.,
∆F = 20 Гц
Пеленгатор: МАСРП
f0 = 2 Гц
Порядок астатизма: 2
Цель: точечная, флюкт.,
∆F = 20 Гц
Порядок выполнения и требования к отчету
В соответствии с номером задания необходимо:
– составить план проведения исследований и согласовать его с
преподавателем,
– произвести необходимые расчеты параметров контура автоматического сопровождения,
– написать собственную программу моделирования в соответствии с заданием и разработанным планом исследований.
В качестве примера используйте файл LAB12.M (см. далее). Рекомендуется создать копию этого файла под своим именем в каталоге \MATLAB\WORK для модификации программы в соответствии
с заданием.
Отчет по лабораторной работе должен включать:
1. Результаты аналитического расчета параметров системы АС и
выбора периода дискретизации Т, коэффициентов передачи цифровых фильтров.
2. Листинг программы и графики, содержащие результаты моделирования.
120
3. Теоретические расчеты зависимостей, исследованных в ходе
выполнения работы.
4. Выводы по работе.
Пример программы моделирования системы
автоматического сопровождения по углу
%
%
%
%
%
%
%
%
%
ОСНОВЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ РЭС
Листинг программы к лабораторной работе #12
Программа моделирования системы автоматического
сопровождения (АС) по углу:
ФНЧ – 1-го порядка,
цель движется с постоянной линейной скоростью V вдоль
окружности радиуса
R, центр окружности расположен в точке стояния РЛС,
ДН парциальных каналов – экспоненциальные.
% Закрываем все открытые окна вывода
close all
% Очищаем рабочую область (WORKSPACE)
clear
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Исходные данные
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Период дискретизации, с
Ts = 1e-3;
N = 1024*8; % Размер выборки
t = (0:N-1)’*Ts; % Текущее время
% Начальное положение цели, град.
alpha_t0 = 0;
% Линейная скорость движения цели, м/с
V = 100;
% Радиус траектории цели, км
R = 100;
% Диаграммы направленности
beta = 2*log(2);
fsum = inline(‘exp(-beta*(alpha/alpha_bw)).^2*cosh(beta*alpha/alpha_
bw)’,’alpha’,’alpha_bw’,’beta’);
fdif = inline(‘exp(-beta*(alpha/alpha_bw)).^2*sinh(beta*alpha/alpha_
bw)’,’alpha’,’alpha_bw’,’beta’);
% Ширина ДН парциального канала, град.
alpha_bw = 1;
% Граничная частота замкнутого контура АС, Гц
F0 = 2;
% Коэффициент перерегулирования
Delta = 0.3;
% Отношение сигнал/шум, дБ
q2 = 20;
121
% Перевод в систему СИ
alpha_t0 = deg2rad(alpha_t0);
R = R*1000;
alpha_bw = deg2rad(alpha_bw);
q2 = 10^(0.1*q2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
Расчет параметров ФАПЧ
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu = beta/alpha_bw; % Крутизна ПХ
xi = 1/sqrt(1 + (pi/log(Delta))^2); % Коэф. переколебательности
W0 = 2*pi*F0; % Перевод в рад/с
tau = 0.5/(xi*W0); % Постоянная времени
k = 0.5*W0/(mu*xi); % Статический коэффициент усиления
% Коэффициенты петлевого фильтра
[numd,dend] = bilinear([0 0 k],[tau 1 0],1/Ts);
% Массив угловых координат цели
alpha_t = rem(alpha_t0 + (V/R)*t,2*pi);
uin = zeros(3,1); % Входной сигнал петлевого фильтра
uout = zeros(2,1); % Выходной сигнал петлевого фильтра
% Массив угловых координат РСН
alpha_a = zeros(N,1);
for n = 2:N
% Суммарный сигнал
esum = fsum(alpha_t(n)-alpha_a(n-1),alpha_bw,beta) + ...
(randn(1) + i*randn(1))/sqrt(0.5*q2); % добавляем
% шум
% Разностный сигнал
edif = fdif(alpha_t(n)-alpha_a(n-1),alpha_bw,beta) + ...
(randn(1) + i*randn(1))/sqrt(0.5*q2); % добавляем
% шум
x = real(edif/esum); % Сигнал ошибки на выходе ФД
uin = [x;uin(1:2)]; % Обновление входного сигнала
% фильтра
y = numd*uin-dend(2:3)*uout; % Выходной сигнал
% фильтра
uout = [y;uout(1)]; % Обновление выходного сигнала
% фильтра
alpha_a(n) = y; % Положение РСН
end
figure(1)
plot(t,alpha_t,t,alpha_a)
grid
xlabel(‘t’)
ylabel(‘\alpha_t(t), \alpha_a(t)’)
legend(‘\alpha_t(t)’, ‘\alpha_a(t)’)
figure(2)
122
plot(t,alpha_t-alpha_a)
grid
xlabel(‘t’)
ylabel(‘\Delta\alpha(t)’)
Контрольные вопросы
1. Каково назначение системы автоматического сопровождения
(АС) по углу?
2. Нарисуйте функциональную схему системы АС по углу и определите назначение основных элементов.
3. Нарисуйте функциональную схему МАСРП и определите назначение основных элементов.
4. Нарисуйте линеаризированную функциональную схему системы АС по углу и определите коэффициенты передачи ее звеньев.
5. Нарисуйте функциональную схема математической модели
системы АС по углу и определите назначение основных элементов.
6. Как рассчитываются параметры математической модели контура автоматического регулирования системы АС по углу?
Литература: [9, с. 76–82; 18]
123
Библиографический список
1. Рабинер, Л. Теория и применение цифровой обработки сигналов / Л. Рабинер, Б. Гоулд; пер. с англ. под ред. Ю. И. Александрова.
М.: Мир, 1978. 848 с.
2. Оппенгейм, А. В. Цифровая обработка сигналов / А. В. Оппенгейм,
Р. В. Шафер; пер. с англ. под ред. С. Я. Шаца. М.: Связь, 1979. 416 с.
3. Антонью, А. Цифровые фильтры: анализ и проектирование /
А. Антонью; пер. с англ. В. А. Лексаченко, В. Г. Челпанова; под ред.
С. А. Понырко. М.: Радио и связь, 1983. 320 с.
4. Сергиенко, А. Б. Цифровая обработка сигналов / А. Б. Сергиенко. СПб.: Питер, 2003. 608 с.
5. Основы цифровой обработки сигналов: курс лекций / А. И. Солонина, Д. А. Улахович, С. М. Арбузов, Е. Б. Соловьева, И. И. Гук.
СПб.: БХВ-Петербург, 2003. 608 с.
6. Монаков, А. А. Основы цифровой обработки сигналов: дискретные сигналы и цифровые фильтры / А. А. Монаков. СПб: ГУАП,
2008. 112 с.
7. Быков, В. В. Цифровое моделирование в статистической радиотехнике / В. В. Быков. М.: Сов. радио, 1971. 328 с.
8. Цифровая обработка сигналов. Моделирование в MATLAB /
А. И. Солонина, С. М. Арбузов. СПб.: БХВ-Петербург, 2008. 816 с.
9. Монаков, А. А. Основы математического моделирования радиотехнических систем / А. А. Монаков. СПб.: ГУАП, 2005. 100 с.
10. Марпл-мл., С. Л. Цифровой спектральный анализ и его приложения / С. Л. Марпл-мл.; пер. с англ. О. И. Хабарова, Г. А. Сидоровой под ред. И. С. Рыжака. М.: Мир, 1990. 584 с.
11. Тихонов, В. И. Оптимальный прием сигналов / В. И. Тихонов.
М.: Радио и связь, 1983. 320 с.
12. Тихонов, В. И. Статистический анализ и синтез радиотехнических систем: учеб. пособие для вузов / В. И. Тихонов, В. Н. Харисов. М.: Радио и связь, 1991. 608 с.
13. Балакришнан, А. Теория фильтрации Калмана / А. Балакришнан; пер. с англ. С. М. Зуева под ред. А. А. Новикова. М.: Мир,
1988. 168 с.
14. Бакалов, В. П. Цифровое моделирование случайных процессов / В. П. Бакалов. М.: САЙНС-ПРЕСС, 2002. 88 с.
15. Кривицкий, Б. Х. Системы автоматической регулировки усиления / Б. Х. Кривицкий, Е. Н. Салтыков. М.: Радио и связь, 1982.
16. Радиоприемные устройства: учеб. пособие для радиотехн.
спец. вузов / Ю. Т. Давыдов, Ю. С. Данич, А. П. Жуковский и др.;
под ред. А. П. Жуковского. М.: Высшая школа, 1989. 342 с.
124
17. Прокис, Дж. Цифровая связь / Дж. Прокис; пер. с англ. Д. Д.
Кловского, Б. И. Николаева; под ред. Д. Д. Кловского. М.: Радио и
связь, 2000. 800 с.
18. Леонов, А. И. Моноимпульсная радиолокация: 2-е изд., перераб. и доп. / А. И. Леонов, К. И. Фомичев. М.: Радио и связь, 1984.
312 с.
125
Содержание
Методические указания................................................ 3
Лабораторная работа № 1. Дискретизация аналоговых
сигналов..................................................................... 4
Лабораторная работа № 2. Исследование линейных
дискретных систем....................................................... 15
Лабораторная работа №3. Автоматизированное
проектирование цифровых фильтров в среде MATLAB...... 28
Лабораторная работа № 4. Методы непараметрического
спектрального анализа случайных стационарных
процессов.................................................................... 41
Лабораторная работа № 5. Параметрический
спектральный анализ случайных стационарных процессов:
метод Юла – Уолкера.................................................... 52
Лабораторная работа № 6. Фильтр Калмана..................... 59
Лабораторная работа № 7. Методы генерации случайных
величин с произвольным законом распределения............. 67
Лабораторные работы № 8 и 9. Моделирование
гауссовских случайных процессов с заданными
спектральными свойствами........................................... 78
Лабораторная работа № 10. Моделирование системы
автоматической регулировки усиления........................... 97
Лабораторная работа № 11. Моделирование системы
фазовой автоматической подстройки частоты.................. 106
Лабораторная работа № 12. Моделирование следящего
моноимпульсного амплитудного суммарно-разностного
пеленгатора................................................................. 113
Библиографический список........................................... 124
126
Документ
Категория
Без категории
Просмотров
4
Размер файла
17 216 Кб
Теги
monakovmirolyubov
1/--страниц
Пожаловаться на содержимое документа