Учебное пособие для студентов высших учебных заведений
 Скачать 5.41 Mb.
|
|
xlabel('Время (с)'); ylabel('X1(t)'); Теперь создадим формирующий фильтр, "пропустим" через него белый шум и результат выведем на рис. 5.37: % Расчет параметров формирующего фильтра 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; % Формирование "профильтрованного" процесса y1=filter(b,a,x1); % Построение графика процесса 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.38 графики модуля ФИ и спектральной плотности мощности (СПМ): Рис. 5.36 5.3. Спектральный и статистический анализ 219 Рис. 5.37 Рис.5.38 % Формирование массива частот df = 1/T; Fmax = 1/Ts; f = - Fmax/2 : df : Fmax/2; dovg = length(f); % Расчет скорректированных массивов Фурье-изображений Fu1 = fft(x1)/dovg; Fu2 = fft(y1)/dovg; Fu1p = fftshift(Fu1); Fu2p = fftshift(Fu2); % Формирование массивов модулей ФИ A1 = abs(Fu1p); A2 = abs(Fu2p); % Вычисление Спектральных Плотностей Мощности S1 = Fu1p.*conj(Fu1p)*dovg; S2 = Fu2p.*conj(Fu2p)*dovg; % Вывод графиков белого шума subplot(2,1,1); stem(f,A1),grid, set(gca,'FontName','Arial Cyr','FontSize',10) title('Модуль ФИ гауссового белого шума'); 5.3. Спектральный и статистический анализ 220 subplot(2,1,2); stem(f,S1),grid, set(gca,'FontName','Arial Cyr','FontSize',10) title('Спектральная плотность мощности'); xlabel('Частота (Гц)'); Рассматривая рис. 5.38, можно убедиться, что спектральная плотность прак- тически одинакова по величине во всем диапазоне частот, чем и обусловлено на- звание процесса - "белый шум". Аналогичную процедуру проведем и с "профильтрованным" процессом (рис. 5.39): Рис. 5.39 Рис. 5.40. % Вывод графиков профильтрованного процесса c1 = fix(dovg/2)-200, c2 = fix(dovg/2)+200, length(f) subplot(2,1,1); stem(f(c1:c2),A2(c1:c2)),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title('Модуль ФИ случайного стационарного процесса'); subplot(2,1,2); stem(f(c1:c2),S2(c1:c2)),grid, set(gca,'FontName','Arial Cyr','FontSize',16) 5.3. Спектральный и статистический анализ 221 title('Спектральная плотность мощности'); xlabel('Частота (Гц)'); Проводя эти вычисления еще раз с новой длительностью процесса T=20 c (см. рис. 5.40), можно наглядно убедиться, что величины ФИ и СПМ практически при этом не изменяются. 5.3.3. Статистический анализ К задачам статистического анализа процессов относятся определение неко- торых статистических характеристик процессов, таких, как корреляционные ха- рактеристики, спектральные плотности мощности и т. д. В предыдущем разделе уже были определены СП случайного процесса на основе установленной связи СП с Фурье-изображением. Однако в Signal Processing Toolbox предусмотрена специальная процедура psd , позволяющая сразу находить СП сигнала. Обращение к ней имеет вид [S,f] = psd (x, nfft, Fmax), где x - вектор заданных значений процесса, - число элементов вектора nfft x , ко- торые обрабатываются процедурой fft , F Ts max = 1 - значение частоты дискрети- зации сигналу, S - вектор значений СП сигнала, f - вектор значений частот, ко- торым соответствуют найденные значения СП. В общем случае длина последних двух векторов равна nfft / 2. Приведем пример применения процедуры psd для нахождения СП преды- дущего случайного процесса: [C,f] = psd(y1,dovg,Fmax); stem(f(1:200),C(1:200)),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title(' Cпектральная плотность мощности'); xlabel('Частота (Гц)'); В результате получим картину, представленную на рис. 5.41 Рис. 5.41 Если ту же процедуру вызвать без указания выходных величин, то результа- том ее выполнения станет выведение графика СП от частоты. 5.3. Спектральный и статистический анализ 222 Рис. 5.42 Например, обращение psd(y1, dovg, Fmax) приведет к построению в графическом окне (фигуре) графика рис. 5.42. При этом значения СП будут откладываться в логарифмическом масштабе в децибелах. Рис. 5.43 5.4. Проектирование фильтров 223 Группа функций xcorr вычисляет оценку Взаимной Корреляционной Функ- ции (ВКФ) двух последовательностей x и y. Обращение c = xcorr(x,y) вычисляет и выдает вектор c длины 2N-1 значений ВКФ векторов x и y длины N. Обращение c = xcorr(x) позволяет вычислить АКФ (автокорреляционную функцию) последо- вательности, заданной в векторе x. Вычислим АКФ для случайного процесса, сформированного ранее: R=xcorr(y1); tau= -10+Ts : Ts : 10; lt=length(tau ); s1r=round(length(R )/2)-lt/2; s2r=round(length(R )/2)+lt/2-1; plot(tau,R(s1r:s2r)),grid set(gca, 'FontName', 'Arial Cyr', 'FontSize', 16) title('АКФ случайного процесса') , xlabel('Запаздывание (с)') На рис. 5.43 представлен результат применения процедуры xcorr. 5.4. Проектирование фильтров 5.4.1. Формы представления фильтров и их преобразования Фильтр как звено системы автоматического управления может быть пред- ставлен в нескольких эквивалентных формах, каждая из которых полностью опи- сывает его: в форме рациональной передаточной функции (tf - представление); если звено является непрерывным (аналоговым), то оно описывается непре- рывной передаточной функцией ) 1 ( ) 2 ( ) 1 ( ) 1 ( ) 2 ( ) 1 ( ) ( ) ( ) ( 1 1 + + + ⋅ + ⋅ + + + ⋅ + ⋅ = = − − n a s a s a m b s b s b s a s b s W n n m m , (5.34) а в случае дискретного фильтра последний может быть представлен дис- кретной передаточной функцией вида n m z n a z a a z m b z b b z a z b z W − − − − ⋅ + + + ⋅ + ⋅ + + + ⋅ + = = ) 1 ( ) 2 ( ) 1 ( ) 1 ( ) 2 ( ) 1 ( ) ( ) ( ) ( 1 1 ; (5.35) в обоих случаях для задания звена достаточно задать два вектора - вектор b коэффициентов числителя и а - знаменателя передаточной функции; в виде разложения передаточной функции на простые дроби; в случае простых корней такое разложение имеет вид (для дискретной передаточ- ной функции): ; ) 1 ( ) 2 ( ) 1 ( ) ( 1 ) ( ) 1 ( 1 ) 1 ( ) ( ) ( ) ( 1 1 1 n m z n m k z k k z n p n r z p r z a z b − − − − − ⋅ + − + + ⋅ + + ⋅ − + + ⋅ − = (5.36) в этой форме звено описывается тремя векторами: вектором-столбцом r вычетов передаточной функции, вектором столбцом p полюсов и векто- ром-строкой k коэффициентов целой части дробно-рациональной функ- ции; 5.4. Проектирование фильтров 224 в каскадной форме (sos - представление), когда передаточная функция звена представлена в виде произведения передаточных функций не выше второго порядка: ∏ ∏ = − − − − = ⋅ + ⋅ + ⋅ + ⋅ + = = L k k k ok k k ok L k k z a z a a z b z b b z H z H 1 2 2 1 1 2 2 1 1 1 ) ( ) ( ; (5.37) параметры каскадного представления задаются в виде матрицы sos, содер- жащей вещественные коэффициенты: ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = L L oL L L oL o o o o a a a b b b a a a b b b a a a b b b sos 2 1 2 1 22 12 2 22 12 2 21 11 1 21 11 1 ; (5.38) в пространстве состояний (ss -представление), т. е. с помощью уравнений звена в форме u D x C y u B x A x ⋅ + ⋅ = ⋅ + ⋅ = & ; (5.39) в этой форме звено задается совокупностью четырех матриц A,B,C и D; путем задания векторов z нулей передаточной функции, p - ее полюсов и k - коэффициента передачи звена (zp - представление): )] ( [ )] 2 ( [ )] 1 ( [ )] ( [ )] 2 ( [ )] 1 ( [ ) ( n p s p s p s m z s z s z s k s W − ⋅⋅ ⋅ − ⋅ − − ⋅⋅ ⋅ − ⋅ − = ; (5.40) решетчатое latc-представление; в этом случае решетчатый фильтр зада- ется векторами k коэффициентов знаменателя решетчатого дискретного фильтра и v - коэффициентов его числителя; коэффициенты k решетчато- го представления некоторого полинома с коэффициентами, представлен- ными вектором а, определяются по этому вектору с помощью рекурсив- ного алгоритма Левинсона. Пакет SIGNAL предоставляет пользователю ряд процедур, позволяющих преобразовать звено (фильтр) из одной формы в другую. Процедуры преобразования к tf-форме 1. Процедура zp2tf осуществляет вычисления векторов b коэффициентов числителя и a знаменателя передаточной функции в форме (5.34) по известным векторам z ее нулей, p - ее полюсов и k - коэффициенту усиления звена. Обраще- ние к процедуре имеет вид: [ b, a] = zp2tf( z, p, k). В общем случае многомерного звена величина z является матрицей, число столбцов которой должно быть равно числу выходов. Вектор-столбец k содер- жит коэффициенты усиления по всем выходам звена. В векторе a выдаются вы- численные коэффициенты знаменателя, а матрица b содержит коэффициенты 5.4. Проектирование фильтров 225 числителей. При этом каждая строка матрицы соответствует коэффициентам чис- лителя для отдельной выходной величины. 2. Процедура ss2tf преобразовывает описание звена (системы) из простран- ства состояний в форму передаточной функции. Обращение к ней вида [ b, a] = ss2tf ( A, B, C, D, iu) позволяет найти коэффициенты числителей (b) и знаменателя (а) передаточных функций системы по всем выходным величинам и по входу с номером "iu", если заданы матрицы A, B, C, D описания системы в виде (5.39). 3. Процедура sos2tf позволяет найти передаточную функцию звена по за- данным параметрам каскадной формы. Для этого надо обратиться к этой процеду- ре таким образом: [ b, a] = sos2tf ( sos), где sos - заданная матрица каскадной формы (5.38). 4. С помощью процедуры latc2tf можно вычислить коэффициенты числите- ля и знаменателя передаточной функции (5.35) по коэффициентам знаменателя и числителя решетчатого фильтра. При этом обращение к ней должно иметь один из видов: [ b, a] = latc2tf( k, v) [ b, a] = latc2tf( k, 'iir') b = latc2tf( k, 'fir') b = latc2tf( k). Первый вид используется, если заданы коэффициенты решетчатого пред- ставления и числителя v и знаменателя k БИХ-фильтра (фильтра с Бесконечной Импульсной Характеристикой). Вторая форма - если решетчатый БИХ-фильтр имеет только полюсы. Третья и четвертая формы применяются для вычисления коэффициентов передаточной функции решетчатого КИХ-фильтра (с Конечной Импульсной Ха- рактеристикой). 5. Нахождение коэффициентов передаточной функции по коэффициентам разложения ее на простые дроби (5.36) осуществляется использованием функций residue и residuez. Первая применяется для непрерывной передаточной функции вида (5.34), вторая - для дискретной передаточной функции (5.35). При обраще- нии вида [ b, a] = residue( r, p, k) [ b, a] = residuez(r, p, k) вычисляются коэффициенты числителя и знаменателя передаточной функции по заданным векторам ее разложения - вычетов r, полюсов p и коэффициентам це- лой части k. С помощью тех же процедур осуществляется разложение заданной пере- даточной функции на простые дроби. При этом обращение к ним должно быть таковым: [ r, p, k] = residue(b, a) [ r, p, k] = residuez(b, a). 5.4. Проектирование фильтров 226 Процедуры перехода от tf - формы к другим 1. Вычисление нулей, полюсов и коэффициентов усиления звена с заданной передаточной функцией можно осуществить, применяя процедуру tf2zp в такой форме: [ z, p, k] = tf2zp(b, a). При этом вектор z будет содержать значения нулей передаточной функции с коэффициентами числителя b и знаменателя а, вектор p - значения полюсов, а k будет равен коэффициенту усиления звена. 2. Нахождение матриц A, B, C и D, описывающих звено с заданной переда- точной функцией в виде совокупности дифференциальных уравнений в форме Коши (5.39), осуществимо с помощью процедуры tf2ss. Если обратиться к ней в форме [A, B, C, D] = tf2ss( b, a), где b и a - соответственно векторы коэффициентов числителя и знаменателя пере- даточной функции, то в результате получим искомые матрицы в указанном по- рядке. 3. Вычисление коэффициентов решетчатого фильтра по заданной дис- кретной передаточной функции можно осуществить при помощи процедуры tf2latc , обращаясь к ней так: [ k, v] = tf2latc( b, a) [ k, v] = tf2latc( 1, a) k = tf2latc( 1, a) k = tf2latc( b). Первое обращение позволяет вычислить коэффициенты k знаменателя и v - числителя решетчатого БИХ-фильтра (модель авторегрессии скользящего средне- го). Обращение во второй форме дает возможность определить вектор коэффици- ентов знаменателя k и скалярный коэффициент v, когда БИХ-фильтр имеет только полюсы (не имеет нулей). В третьей форме определяются только коэффициенты знаменателя решетчатого фильтра. Наконец, четвертая форма предназначена для нахождения вектора k коэффициентов решетчатого КИХ-фильтра (задаваемого только вектором b коэффициентов числителя передаточной функции). Другие преобразования 1. Вычисление коэффициентов решетчатого представления по коэффици- ентам полинома можно осуществить, используя функцию poly2rc. Обращение k = poly2rc( а) позволяет найти коэффициенты решетчатого представления k по коэффициентам а заданного полинома. Вектор а должен содержать только вещественные элемен- ты и должно выполняться условие 0 ≠ a . Размер вектора k на единицу меньше размера вектора а коэффициентов полинома. По коэффициентам решетчатого представления очень просто определить, находятся ли все полюсы внутри единичного круга. Для этого достаточно прове- 5.4. Проектирование фильтров 227 рить, что все элементы вектора k по абсолютной величине не превышают едини- цы. Обратная задача вычисления коэффициентов полинома по коэффициентам решетчатого представления решается путем применения функции rc2poly: a = rc2poly(k). 2. Процедура sos2ss определяет матрицы (5.39) A, B. C и D, описывающие звено в пространстве состояний, по заданной матрице SOS каскадной формы (5.38): [A,B,C,D] = sos2ss (SOS). Элементы матрицы SOS должны быть вещественными. Обратный переход осуществляется при помощи функции ss2sos SOS = ss2sos(A,B,C,D). 3. Функция sos2zp дает возможность определить (5.40) векторы z нулей, p - полюсов и коэффициент усиления k звена, заданного каскадной формой переда- точной функции (т. е. матрицей SOS (5.38)): [z,p,k] = |