Учебное пособие для студентов высших учебных заведений
 Скачать 5.41 Mb.
|
|
plot(t,y), grid, set(gca,'FontName','Arial Cyr','FontSize',16) title(' Пример процедуры CHIRP') xlabel('Время ( с )'), ylabel('Выходной процесс Y(t)') Результат показан на рис. 5.12. Рис. 5.12. Еще одна процедура diric формирует массив значений так называемой функции Дирихле, определяемой соотношениями: diric t t k k nt n t п и д угих t k n ( ) , , , sin( / ) sin( / ) р р ( ) = − = = ± ⋅ ⎧ ⎨ ⎪ ⎩ ⎪ − 1 2 0 1 2 2 1 π ,.. ± 2 Функция Дирихле является периодической. При нечетных n период равен 2 π , при четных - 4 π . Максимальное значение ее равно 1, минимальное -1. Параметр n должен быть целым положительным числом. Обращение к функции имеет вид: y = diric(t, n). Ниже приведены операторы, которые иллюстрируют использование процедуры diric и выводят графики функции Дирихле для n =3, 4 и 5: t=0 : 0.01 : 50; y1=0.7*diric(pi*t/5, 3); 5.1. Формирование типовых процессов 198 plot(t,y1), grid,set(gca,'FontName','Arial Cyr','FontSize',16) title('Функция Дирихле Y(t) = 0.7* DIRIC(pi*t/5, 3)') xlabel('Время ( с )'), ylabel('Выходной процесс Y(t)') Результаты представлены на рис. 5.13...5.15. Рис. 5.13. Рис. 5.14. Рис. 5. 15. 5.2. Общие средства фильтрации 199 5.2. Общие средства фильтрации. Формирование слу- чайных процессов 5.2.1. Основы линейной фильтрации Рассмотрим основы линейной фильтрации на примере линейного стацио- нарного фильтра, который в непрерывном времени описывается дифференциаль- ным уравнением второго порядка: && & y y y o o + ⋅ + ⋅ = 2 2 ζω ω A x ⋅ , (5.1) где x - заданный процесс, подаваемый на вход этого фильтра второго порядка; - процесс, получаемый на выходе фильтра; y ω o - частота собственных колебаний фильтра, а ζ - относительный коэффициент затухания этого фильтра. Передаточная функция фильтра, очевидно, имеет вид: W s y s x s A s s o o ( ) ( ) ( ) = = + ⋅ + 2 2 ζω ω 2 ) (5.2) Для контроля и графического представления передаточной функции любого линейного динамического звена удобно использовать процедуру freqs. В общем случае обращение к ней имеет вид: h = freqs(b, a, w). При этом процедура создает вектор h комплексных значений частотной ха- рактеристики W j ( ω по передаточной функции W s ( ) звена, заданной векторами коэффициентов ее числителя b и знаменателя a, а также по заданному век-тору w частоты ω . Если аргумент w не указан, процедура автоматически выбирает 200 отсчетов частоты, для которых вычисляется частотная характеристика. Если не указана выходная величина, т.е. обращение имеет вид freqs(b, a, w), процедура выводит в текущее графическое окно два графика - АЧХ и ФЧХ. Приведем пример. Пусть для передаточной функции (5.2) выбраны такие значения параметров: A T o o = = = = 1 0 05 2 ; . ; / 1 ζ π ω Вычислим значения коэффициентов числителя и знаменателя и выведем графики АЧХ и ФЧХ: T0=1; dz=0.05; om0=2*pi/T0; A=1; a1(1)=1; a1(2)=2*dz*om0; a1(3)= om0^2; b1(1)=A; freqs(b1,a1) В результате получим рис. 5.16. Допустим, что заданный процесс x (t) представлен в виде отдельных его значений в дискретные моменты времени, которые разделены одинаковыми про- межутками Ts времени (дискретом времени). Обозначим через x k ( ) значение процесса в момент времени t k T = s ⋅ , где - номер измерения с начала процесса. k 5.2. Общие средства фильтрации 200 Рис. 5.16 Запишем уравнение (5.1) через конечные разности процессов x и , учи- тывая, что конечно-разностным эквивалентом производной является конечная разность y &y ∆y k Ts y k y k Ts ( ) ( ) ( ) = − −1 , а эквивалентом производной второго порядка является конечная разность вто- рого порядка &&y ∆ ∆ ∆ 2 2 2 2 1 2 1 y k Ts y k y k Ts y k y k y k Ts ( ) ( ) ( ) ( ) ( ) ( ) = − − = − − + − 2 ( ) Тогда разностное уравнение ( ) ( ) ( ) ( ) ( ) 1 2 2 1 1 2 2 2 2 + ⋅ + ⋅ ⋅ − + ⋅ ⋅ − + − = ⋅ ⋅ ζω ω ζω o o o Ts Ts y k Ts y k y k A Ts x k (5.3) является дискретным аналогом дифференциального уравнения (5.1). Применяя к полученному уравнению Z-преобразование, получим: y z a a z a z A T s x z o ( ) [ ] ( ) ⋅ + ⋅ + ⋅ = ⋅ ⋅ − − 1 1 2 2 2 , (5.4) где a T s o o o = + ⋅ + ⋅ 1 2 2 2 ζω ω ; T s Ts a o 1 2 1 = − + ⋅ ( ); ζω (5.5) a 2 1 = Дискретная передаточная функция фильтра определяется из уравнения (5.4): 5.2. Общие средства фильтрации 201 G z y z x z A T s a a z a z o ( ) ( ) ( ) = = ⋅ + ⋅ + ⋅ − 2 1 1 2 2 − ) , (5.6) Таким образом, цифровым аналогом ранее введенного колебательного звена является цифровой фильтр с коэффициентами числителя и знаменателя, рассчи- танными по формулам (5.4) и (5.5): T0=1; dz=0.05; Ts=0.01; om0=2*pi/T0; A=1; oms=om0*Ts; a(1)= 1+2*dz*oms+oms^2; a(2)= - 2*(1+dz*oms); a(3)=1; b(1)=A*Ts*Ts*(2*dz*om0^2); Чтобы получить частотную дискретную характеристику G e j ( ω по дис- кретной передаточной функции G z , которая задана векторами значений ее чис- лителя b и знаменателя a, удобно использовать процедуру freqz, обращение к ко- торой аналогично обращению к процедуре freqs : ( ) freqz( b,a) Рис. 5.17 В системе MatLAB фильтрация, т.е. преобразование заданного сигнала с помощью линейного фильтра, описываемого дискретной передаточной функцией вида G z y z x z b b z b z a a z a z o m o n n ( ) ( ) ( ) = = + ⋅ + + ⋅ + ⋅ + + ⋅ m − − − 1 1 1 1 − (5.7) осуществляется процедурой filter следующим образом y = filter(b, a, x) , где x - заданный вектор значений входного сигнала; y - вектор значений выходно- го сигнала фильтра, получаемого вследствие фильтрации; b - вектор коэффициен- 5.2. Общие средства фильтрации 202 тов числителя дискретной передаточной функции (5.7) фильтра; a - вектор коэф- фициентов знаменателя этой функции. В качестве примера рассмотрим такую задачу. Пусть требуется получить достаточно верную информацию о некотором "полезном" сигнале, имеющем си- нусоидальную форму с известным периодом Т 1 =1с и амплитудой А 1 =0.75. Сформируем этот сигнал как вектор его значений в дискретные моменты времени с дискретом Тs = 0.001 с: Ts=0.001; t=0 : Ts : 20; A1=0.75; T1=1; Yp=A1*sin(2*pi*t/T1); plot(t(10002:end),Yp(10002:end)),grid,set( gca,'FontName','ArialCyr','FontSize',16) title('"Полезный" процесс '); xlabel('Время (с)'); ylabel('Yp(t)') Рис. 5.18 Рис. 5.19 5.2. Общие средства фильтрации 203 Допустим, что вследствие прохождения через ПП (первичный преобразова- тель) к полезному сигналу добавился шум ПП в виде более высокочастотной си- нусоиды с периодом Т 2 = 0.2с и амплитудой А 2 =5, а в результате измерения к не- му еще добавился белый гауссовый шум измерителя с интенсивностью А ш =5. В результате создался такой измеренный сигнал x t ( ) (см. рис. 5.19): T2=0.2;; A2=10; eps=pi/4; Ash=5; x=A1.*sin(2*pi*t./T1)+A2.*sin(2*pi.*t./T2+eps)+Ash*randn(1,length(t)); plot(t(10002:end),x(10002:end)),grid,set(gca,'FontName','ArialCyr','FontSize',16), title('Входной процесс '); xlabel('Время (с)'); ylabel('X(t)') Требуется так обработать измеренные данные x, чтобы восстановить по ним полезный процесс как можно точнее. Так как частота полезного сигнала заранее известна, восстановление его можно осуществить при помощи резонансного фильтра отмеченного выше вида. При этом необходимо создать такой фильтр, чтобы период его собственных коле- баний Т ф был равен периоду колебаний полезного сигнала (Т ф = Т 1 ). Для того, чтобы после прохождения через такой фильтр амплитуда восстановленного сиг- нала совпадала с амплитудой полезного сигнала, нужно входной сигнал фильтра умножить на постоянную величину 2 (ибо при резонансе амплитуда выходно- го сигнала "уменьшается" именно во столько раз по сравнению с амплитудой входного сигнала). Сформируем фильтр, описанный выше: 2 ςω o T1=1; Tf=T1; dz=0.05; om0=2*pi/Tf; A=1; oms=om0*Ts; a(1)= 1+2*dz*oms+oms^2; a(2)= - 2*(1+dz*oms); a(3)=1; b(1)=A*Ts*Ts*(2*dz*om0^2); и "пропустим" сформированный процесс через него y=filter(b,a,x) plot(t(10002:end),y(10002:end),t(10002:end),Yp(10002:end)),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title('Процесс на выходе фильтра (Tf=1; dz=0.05)'); xlabel('Время (с)'); ylabel('Y(t)') Рис. 5.20 5.2. Общие средства фильтрации 204 В результате получаем восстановленный процесс (рис. 5.20). Для сравнения на этом же графике изображен восстанавливаемый процесс. Как видим, созданный фильтр достаточно хорошо восстанавливает полез- ный сигнал. Однако более точному восстановлению препятствуют два обстоя- тельства: 1) восстановленный процесс устанавливается на выходе фильтра только спустя некоторое время вследствие нулевых начальных условий самого фильтра как динамического звена; это иллюстрируется ниже, на рис. 5.21; y=filter(b,a,x) plot(t,y,t,Yp),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title('Процесс на выходе фильтра (Tf=1; dz=0.05)'); xlabel('Время (с)'); ylabel('Y(t)') Рис. 5.21 2) в установившемся режиме наблюдается значительный сдвиг ( π / 2 ) фаз между восстанавливаемым и восстановленным процессами; это тоже понятно, так как при резонансе сдвиг фаз между входным и выходным процессами достигает именно такой величины. Чтобы избежать фазовых искажений полезного сигнала при его восстанов- лении, можно воспользоваться процедурой двойной фильтрации - filtfilt. Обраще- ние к ней имеет такую же форму, что и к процедуре filter. В отличие от последней, процедура filtfilt осуществляет обработку вектора x в два приема: сначала в пря- мом, а затем в обратном направлении. Результат применения этой процедуры в рассматриваемом случае приве- ден на рис. 5.22. y=filtfilt(b,a,x) plot(t,y,t,Yp),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title('Применение процедуры FILTFILT'); xlabel('Время (с)'); ylabel('Y(t)') 5.2. Общие средства фильтрации 205 Рис. 5.22 5.2.2. Формирование случайных процессов В соответствии с теорией, сформировать случайный процесс с заданной корреляционной функцией можно, если сначала сформировать случайный про- цесс, являющийся нормально (по гауссовому закону) распределенным белым шу- мом, а затем "пропустить" его через некоторое динамическое звено (формирую- щий фильтр). На выходе получается нормально распределенный случайный про- цесс с корреляционной функцией, вид которой определяется типом формирующе- го фильтра как динамического звена. Рис. 5.23 Белый гауссовый шум в MatLAB образуется при помощи процедуры randn. Для этого достаточно задать дискрет времени Ts, образовать с этим шагом массив 5.2. Общие средства фильтрации 206 (вектор) t моментов времени в нужном диапазоне, а затем сформировать по ука- занной процедуре вектор-столбец длиною, равной длине вектора t, например Ts=0.01; t=0 : Ts : 20; x1=randn(1,length(t)); Построим график полученного процесса: plot(t,x1),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title('Входной процесс - белый шум Гаусса (Ts=0.01)'); xlabel('Время (с)'); ylabel('X1(t)') Процесс x 1 (t) с дискретом времени Ts=0.01c представлен на рис. 5.23. Для другого значения дискрета времени (Ts=0.001c), повторяя аналогичные операции, получим процесс x 2 (t), изображенный на рис. 5.24. Рис. 5.24 Создадим дискретный фильтр второго порядка с частотой собственных ко- лебаний ω π o = 2 рад./с = 1 Гц и относительным коэффициентом затухания ζ = 0 05 по формулам (5.5) коэффициентов: om0=2*pi; dz=0.05; A=1; oms=om0*Ts; a(1)= 1+2*dz*oms+oms^2; a(2)= - 2*(1+dz*oms); a(3)=1; b(1)=A*2*dz*oms^2; "пропустим" образованный процесс x 1 (t) через созданный фильтр: y1=filter(b,a,x1); и построим график процесса y 1 (t) на выходе фильтра: plot(t,y1),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title('Процесс на выходе фильтра (T0=1; dz=0.05, Ts= 0.01)'); xlabel('Время (с)'); ylabel('Y1(t)') Результат представлен на рис. 5.25. Аналогичные операции произведем с процессом x 2 (t). В результате полу- чим процесс y 2 (t), приведенный на рис. 5.26. Ts=0.001; om0=2*pi; dz=0.05; A=1; oms=om0*Ts; a(1)= 1+2*dz*oms+oms^2; a(2)= - 2*(1+dz*oms); a(3)=1; b(1)=A*2*dz*oms^2; 5.2. Общие средства фильтрации 207 y1=filter(b,a,x2); t=0 : Ts : 20; plot(t,y1),grid, set(gca,'FontName','Arial','FontSize',14) title('Процес на виході фільтра (T0=1; dz=0.05, Ts= 0.001)'); xlabel('Час (с)'); ylabel('Y2(t)') Рис. 5.25. Рис. 5.26. Как видим, на выходе формирующего фильтра действительно образуется случайный колебательный процесс с преобладающей частотой 1 Гц. 5.3. Спектральный и статистический анализ 208 5.3. Процедуры спектрального (частотного) и статистического анализа процессов 5.3.1. Основы спектрального и статистического анализа Основная задача спектрального анализа сигналов - выявление гармониче- ского спектра этих сигналов, т.е. определение частот гармонических составляю- щих сигнала (выявление частотного спектра), амплитуд этих гармонических со- ставляющих (амплитудного спектра) и их начальных фаз (фазового спектра). В основе спектрального анализа лежит теория Фурье о возможности разло- жения любого периодического процесса с периодом T f = = 2 1 π ω (где ω - круго- вая частота периодического процесса, а f - его частота в герцах) в бесконечную, но счетную сумму отдельных гармонических составляющих. Напомним некоторые положения спектрального анализа. Прежде всего, любой периодический процесс с периодом T может быть представлен в виде так называемого комплексного ряда Фурье: x t X m e X m e j mf t m j m t m ( ) ( ) ( ) * ( ) * ( = ⋅ = ⋅ ⋅ =−∞ +∞ ⋅ =−∞ +∞ ∑ ∑ 2 π ) ω , (5.8) причем комплексные числа , которые называют X m * ( ) комплексными амплиту- дами гармонических составляющих, вычисляются по формулам: X m T x t e dt T x t e dt j mf t T T j m t T T * ( ) ( ) ( ) ( ) = ⋅ = ⋅ − − − − ∫ ∫ 1 1 2 2 2 2 2 π ( ) ω . (5.9) Таким образом, частотный спектр периодического колебания состоит из частот, кратных основной (базовой) частоте f , т.е. частот f m f m m = ⋅ = ( , , ,. 0 1 2 ..) (5.10) Действительные и мнимые части комплексных амплитуд образуют соответственно X m * ( ) действительный и мнимый спектры периодического колебания. Если комплексную амплитуду (5.9) представить в экспоненциальной форме X m a e m j m * ( ) = ⋅ ⋅ 2 ϕ , (5.11) то величина будет представлять собой амплитуду гармонической составляю- щей с частотой a m f m m f = ⋅ , а ϕ m - начальную фазу этой гармоники, имеющей форму косинусоиды, т. е. исходный процесс можно также записать в виде x t a a mft o m m ( ) cos( ) = + ⋅ + = +∞ ∑ 2 1 π ϕ m , (5.12) который, собственно, и называют рядом Фурье. Для действительных процессов справедливы следующие соотношения |