Учебное пособие для студентов высших учебных заведений
 Скачать 5.41 Mb.
|
|
1.4.5. Векторная фильтрация и спектральный анализ В системе MatLAB есть несколько функций для проведения цифрового ана- лиза данных наблюдений (измерений). Так, функция y = filter(b,a,x) обеспечивает формирование вектора y по за- данным векторам b, a, x в соответствии с соотношением: y(k) = b(1)*x(k) + b(2)*x(k-1) + ... + b(nb+1)*x(k-nb) - -a(2)*y(k-1) - a(3)*y(k-3) - ... - a(na+1)*y(k-na), (1.1) где вектор b имеет такой состав b = [ b(1), b(2),... , b(nb+1)], a вектор а a = [1, a(2), a(3),... , a(na+1)]. Соотношение (1) можно рассматривать как конечно-разностное уравнение фильтра с дискретной передаточной функцией вида рациональной дроби, коэф- фициенты числителя которого образовывают вектор b, а знаменателя - вектор а, на вход которого подается сигнал x(t), а на выходе формируется сигнал y(t). Тогда вектор y будет представлять собой значение исходного сигнала этого фильтра в дискретные моменты времени, соответствующие заданным значениям входного сигнала x(t) (вектор х). Ниже приведен пример применения функции filter. » x = 0:0.1:1; » b = [1 2]; » a = [ 1 0.1 4]; » y = filter(b,a,x) y = Columns 1 through 7 0 0.1000 0.3900 0.2610 -0.5861 0.3146 3.9129 Columns 8 through 11 0.2503 -13.4768 2.8466 56.4225 1.4. Функции прикладной численной математики 67 Функции fft (Fast Fourier Transformation) и ifft (Invers Fast Fourier Transformation) осуществляют преобразование заданного вектора, сооствествую- щие дискретному прямому и обратному преобразованиям Фурье. Обращение к этим функциям вида: y = fft ( x, n ); x = ifft ( y, n ) приводит к формированию вектора y в первом случае и х - в втором по форму- лам: ∑ = − ⋅ − ⋅ ⋅ − ⋅ = n m n k m j e m x k y 1 / ) 1 ( ) 1 ( 2 ) ( ) ( π ; (1.2) ∑ = − ⋅ − ⋅ ⋅ ⋅ = n k n k m j e k y n m x 1 / ) 1 ( ) 1 ( 2 ) ( 1 ) ( π , (1.3) где j - обозначение мнимой единицы; n - число элементов заданного вектора х (оно представляет также размер выходного вектора y). Приведем пример. Сформируем входной сигнал в виде вектора, элементы которого равняются значениям функции, являющейся суммой двух синусоид с частотами 5 и 12 Гц (рис. 1.21). t =0:0.001:2; x = sin(2*pi*5*t) + cos(2*pi*12*t); plot(t, x); grid set(gcf,'color','white') set(gca,'FontName','Arial Cyr','FontSize',16), title('Входной процесс '); xlabel('Время (с)'); ylabel('Х(t)') 0 0.5 1 1.5 2 -2 -1 0 1 2 Входной процесс Время (с) Х (t ) Рис. 1.21 Найдем Фурье-изображение этого сигнала и выведем графическое пред- ставление модуля його Фурье-изображения: y = fft(x); 1.4. Функции прикладной численной математики 68 a =abs(y); plot(a); grid set(gca,'FontName','Arial Cyr','FontSize',16), title('Модуль Фурьє - изображения '); xlabel('Номер элемента вектора'); ylabel('abs(F(X(t))') Результат отображен на рис. 1.22. 0 500 1000 1500 2000 2500 0 200 400 600 800 1000 1200 Модуль Фурьє - изображения Номер элемента вектора a bs (F (X (t )) Рис. 1.22 Теперь осуществим обратное преобразование с помощью функции ifft и результат также выведем в форме графика: z = ifft(y); plot(t, z); grid set(gca,'FontName','Arial Cyr','FontSize',16), title('Обратное Фурье-преобразование '); xlabel('Время (с)'); ylabel('Z(t)') На рис. 1.23 изображен результат. Рассматривая его, можно убедиться, что воспроизведенный процесс точно совпадает с исходным. Внимательно изучая формулу дискретного преобразования Фурье, можно заметить: а) номер m отвечает моменту времени t m , в который измерен входной сиг- нал x(m); при этом ; t 1 0 = б) номер k - это индекс значения частоты f k , которому отвечает найденный элемент y(k)дискретного преобразования Фурье; в) чтобы перейти от индексов к временной и частотной областям, необхо- димо знать значение h дискрета (шага) времени, через который измерен входной сигнал x(t) и промежуток T времени, на протяжении которого он измеряется; то- гда шаг (дискрет) по частоте в изображении Фурье определится соотношением: Df = 1/T, (1.4) а диапазон изменения частоты - формулой 1.4. Функции прикладной численной математики 69 F = 1/h; (1.5) так, в анализируемом примере (h =0. 001, T = 2, n = 21): Df = 0.5; F = 1000; 0 0.5 1 1.5 2 -2 -1 0 1 2 3 Обратное Фурье-преобразование Время (с) Z( t) Рис. 1.23 г) из (2) следует, что индексу k = 1 отвечает нулевое значение частоты(f 0 = 0); иначе говоря, первый элемент вектора y(1) является значением Фурье- изображения при нулевой частоте, т. е. - просто суммой всех заданных значений вектора x; отсюда получаем, что вектор y(k) содержит значение Фурье- изображения, начиная из частоты f 0 = 0 (которой отвечает k = 1) до максимальной частоты fmax = F (которой отвечает k = n); таким образом, Фурье-изображение определяется функцией fft только для положительных частот в диапазоне от 0 к F; это неудобно для построения графиков Фурье-изображения от частоты; более удобным и привычным представляется переход к вектору Фурье-изображения, определенному в диапазоне частот от (-F/2) до F/2; частота получила название F F N = / 2 частоты Найквиста; д) как известно, функция e j z ⋅ является периодической по z с периодом 2 π ; поэтому информация об Фурье-изображении при отрицательных частотах расположена во второй половине вектора y(k). Сформируем для анализируемого примера массив частот, исходя из выше- упомянутого: f = 0 : 0.5 : 1000; и выведем графік с аргументом-частотой (рис. 1.24): plot(f,a); grid set(gca,'FontName','Arial Cyr','FontSize',16), title('Модуль Фурье - изображения '); xlabel('Частота (Гц)'); ylabel('abs(F(X(t))') 1.4. Функции прикладной численной математики 70 Как следует из рассмотрения рис. 1.24, по нему непросто распознать те час- тоты (5 и 12 Гц), с которыми изменяется входной сигнал. Это - следствие того об- стоятельства, которое было отмечено в примечании г). Чтобы определить частот- ный спектр входного сигнала, нужно сначала преобразовать полученный вектор y Фурье-изображения с помощью процедуры fftshift. 0 200 400 600 800 1000 0 200 400 600 800 1000 1200 Модуль Фурье - изображения Частота (Гц) a bs (F (X (t )) Рис. 1.24 Функция fftshift (обращение к ней осуществляется таким образом: z = fftshift (y)) предназначена для формирования нового вектора z из заданного век- тора y путем перестановки второй половины вектора y в первую половину век- тора z . При этом вторая половина вектора z состоит из элементов первой полови- ны вектора y. Более точно эту операцию можно задать соотношениями: z(1) = y(n/2+1); ... , z(k) = y(n/2+k);... , z(n/2) = y(n); z(n/2+1) = y(1);... ..., z(n/2+k) = y(k);... z(n) = y(n/2). Примечание. Операцию fftshift удобно использовать для преобразования массива Фурье-изображение с целью построения его графика в частотной об- ласти. Тем не менее этот массив не может быть использован для обратного преобразования Фурье. Проиллюстрируем применение этой функции к предыдущему примеру: f1 = -500 : 0.5 : 500; % Перестройка вектора частот v = fftshift(y); % Перестройка вектора Фурье-изображения a =abs(v); % Отыскание модуля % Вывод графика plot(f1(970:1030),a(970:1030)); grid set(gca,'FontName','Arial Cyr','FontSize',16), title('Модуль Фурье - изображения'); xlabel('Частота (Гц)'); ylabel('abs(F(X(t))') 1.4. Функции прикладной численной математики 71 -20 -10 0 10 20 0 200 400 600 800 1000 1200 Модуль Фурье - изображения Частота (Гц) a bs (F (X (t )) Рис. 1.25 Из графика рис. 1.25 уже становится очевидным, что в спектре входного сигнала есть две гармоники - с частотами 5 и 12 Гц. Остается лишь то неудобство, что из графика спектра невозможно устано- вить амплитуды этих гармоник. Во избежание этого, нужно весь вектор y Фурье- изображения разделить на число его элементов (n), чтобы получить вектор ком- плексного спектра сигнала: N=length(y); a =abs(v)/N; plot(f1(970:1030),a(970:1030)); grid set(gca,'FontName','Arial Cyr','FontSize',16,'Color','white'), title('Модуль комплексного спектра'); xlabel('Частота (Гц)'); ylabel('abs(F(X(t)) / N') Результат приведен на рис. 1.26. Из его рассмотрения вытекает, что "ам- плитуды" всех составляющих гармоник равны 0.5. При этом нужно принять во внимание, что "амплитуды" распределены между положительными и отрицатель- ными частотами поровну, поэтому они вдвое меньшее действительной амплитуды соответствующей гармоники. 1.4.6. Задания Задание 1.6. 1. Введите произвольную матрицу размером (4*6). Найдите сумму наи- больших элементов ее строк. 2. Введите квадратную матрицу (5*5) с одним наименьшим элементом. Найдите сумму элементов строки, в которой размещен элемент с наименьшим значением. 1.4. Функции прикладной численной математики 72 -20 -10 0 10 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Модуль комплексного спектра Частота (Гц) a bs (F (X (t )) / N Рис. 1.26 3. Введите матрицу (6*9), в которой есть единственные наибольший и наи- меньшие элементы и они расположены в разных строках. Поменяйте местами строку с наибольшим элементом и строку с наименьшим элементом. 4. Введите матрицу (5*6) с разными значениями элементов. В каждой стро- ке выберите элемент с наименьшим значением, из полученных чисел выберите наибольшее. Найдите индексы полученных элементов. 5. Введите матрицу (5*6). Найдите вектор, элементами которого являются наибольшие элементы соответствующей строки матрицы. 6. Введите матрицу (5*6). Постройте вектор, элементами которого являются суммы наибольшего и наименьшего элементов соответствующей строки матрицы. 7. Введите матрицу (5*6). Постройте вектор, элементами которого являются средние значения элементов соответствующей строки матрицы. 8. Введите матрицу (5*6). Постройте вектор, элементами которого являются среднеквадратичные отклонения элементов соответствующей строки матрицы от их среднего значения. 9. Введите матрицу (5*6). Постройте вектор, элементами которого являются средние арифметические наибольшего и наименьшего элементов соответствую- щей строки матрицы. 10. Введите матрицу (6*5). Постройте вектор, элементами которого являют- ся суммы квадратов элементов соответствующего столбца матрицы. 11. Введите матрицу (5*5). Постройте векторы, элементами которых явля- ются суммы элементов столбцов матрицы, произведения элементов столбцов и наименьшие значения элементов столбцов. 12. Введите матрицу (5*6). Найдите среднее арифметическое наибольших и наименьших ее элементов. 13. Введите матрицу (5*5). Постройте вектор, элементами которого являют- ся элементы главной диагонали матрицы. Найдите след матрицы. 1.4. Функции прикладной численной математики 73 14. Введите две матрицы (4*4). Постройте новую матрицу размером (4*8), включая в первые 4 столбца строки первой матрицы, а в другие - столбцы второй матрицы. 15. Найдите сумму всех элементов матрицы размером (4*3). Задание 1.7. Вычислите векторы: а) модуля частотной передаточной функции (ЧПФ); б) аргумента ЧПФ; в) действительной части ЧПФ; г) мнимой части ЧПФ по заданным числителю и знаменателю передаточной функции (таблица 1.4). Предварительно найдите корни знаменателя передаточной функции, опре- делите наибольшую собственную частоту ω max системы. Обеспечьте вычисление ЧПФ при 100 значениях частоты ω в диапазоне от 0 до 5ω max Таблица 1.4 Ва- ри- ант Числитель Знаменатель 1 1. 82p+67.56 p 4 +2. 65p 3 +3. 09p 2 +7. 04p+34.05 2 4. 61p 2 +1. 82p+67.56 p 4 +3. 65p 3 +45p 2 +7. 04p+125 3 p 2 +4p+23 p 4 +2p 3 +39p 2 +2p+45 4 3p 2 +1. 82p+67.56 p 2 +7. 04p+34.05 5 p+6 p 2 +0. 7p+48 6 p 3 +4. 61p 2 +1. 82p 2. 65p 3 +3p 2 +4p+87 7 p 3 +4. 61p 2 +1. 82p+67.56 p 4 +2. 65p 3 +68p 2 +5p+34 8 4. 61p 2 +68 p 4 +2. 65p 3 +3. 09p 2 +7. 04p+34.05 9 7.56 p 4 +2. 65p 3 +3. 09p 2 +7. 04p+34.05 10 p 3 +1. 8p+7 p 4 +6. 5p 3 +39p 2 +7p+45 11 p 3 +4. 61p 2 +1. 82p+67.56 p 3 +3. 09p 2 +70p+34 12 p 2 +1. 8p+78 2. 65p 3 +3. 09p 2 +7. 04p+34.05 13 p 3 +1. 82p+67.56 p 4 +2. 6p 3 +3p 2 +4p+34 14 p 3 +4. 61p 2 +1. 82p+67.56 7p 2 +7p+34 15 4. 61p 2 +1. 82p+67.56 p 2 +7. 04p+560 16 1. 82p+67.56 3. 09p 2 +7. 04p+34.05 17 p 3 3. 09p 2 +7. 8p+125 18 1. 82p p 3 +3. 09p 2 +7. 04p+34.05 19 4. 61p 2 p 2 +7. 04p+34.05 20 p 3 +67.56 p 4 +2. 65p 3 +3. 09p 2 +7. 04p+34.05 21 p 3 p 4 +2p 3 +3p 2 +12p+100 22 p 3 +4. 61p 2+ 1. 82p+67.56 p 4 +5p 3 +30p 2 +7p+305 23 p 2 +1. 82p+67.56 p 4 +2p 3 +9p 2 +4p+35 24 p 3+ 61p 2 +182p+67 p 4 +3p 3 +9p 2 +0. 04p+39 1.4. Функции прикладной численной математики 74 25 p 2 +1. 82p+67.56 p 4 +5p 3 +20p 2 +7p+34 Указание. Частотной передаточной функцией называют передаточную функцию системы при мнимых значениях її аргумента ( p j = ⋅ ω ). Собственные частоты системы - это значения модулей мнимых частей корней характеристического уравнения системы (которое получается прирав- ниванием нулю знаменателя передаточной функции). Задача 1.8. Введите произвольную матрицу размером (5*5). Найдите: 1) определитель матрицы; в случае, если определитель окажется равным нулю, или слишком малым, измените некоторые элементы матрицы и по- вторите вычисления; 2) обратную матрицу; проверьте правильность путем обращения обратной матрицы; 3) характеристический полином матрицы; 4) корни характеристического полинома матрицы; рассортируйте корни по комплексно-спряженным парам и в порядке возрастания величин; 5) собственные значения матрицы; сравните с ранее найденными корнями характеристического полинома; 6) LU-разложение матрицы; проверьте его правильность; 7) QR-разложение матрицы; проверьте его правильность; 8) сингулярные числа матрицы; сравните их с получаемыми при svd- разложении; 9) след матрицы; 10) число обусловленности матрицы; 11) экспоненту от матрицы; 12) логарифм от экспоненты матрицы; сравните с исходной матрицей. 1.4.7. Вопросы 1. Какой объект в MatLAB называется полиномом? 2. Как в MatLAB осуществляется перемножение и деление полиномов? 3. При помощи каких функций можно найти корни заданного полинома, значение полинома по известному значению аргумента? 4. Какие функции позволяют найти производную от полинома? 5. Как найти характеристический полином матрицы? 1.5. Построение простейших графиков 75 1.5. Построение простейших графиков 1.5.1. Процедура plot Вывод графиков в системе MatLAB - настолько простая и удобная опера- ция, что ее можно использовать даже в режиме калькулятора. Основной функцией, обеспечивающей построение графиков на экране дис- плея, является функция plot. Общая форма обращения к ней такова: plot(x1, y1, s1, x2, y2, s2,...). Здесь x1, y1 - заданные векторы, элементами которых являются массивы значений аргумента (х1) и функции (y1), отвечающие первой кривой графика; x2, y2 - массивы значений аргумента и функции второй кривой и т.д. При этом пред- полагается, что значения аргумента откладываются вдоль горизонтальной оси графика, а значения функции - вдоль вертикальной оси. Переменные s1, s2,... яв- ляются символьными (их указание не является обязательным). Любая из них мо- жет содержать до трех специальных символов, определяющих соответственно: а) тип линии, которая соединяет отдельные точки графика; б) тип точки графика; в) цвет линии. Если переменные s не указаны, то тип линии по умолчанию - отре- зок прямой, тип точки - пиксел, а цвет устанавливается (в версии 5.3) в такой очередности: - синий, зеленый, красный, голубой, фиолетовый, желтый, черный и белый - в зависимости от того, какая по очереди линия выводится на график. Например, обращение вида |