Учебное пособие для студентов высших учебных заведений
 Скачать 5.41 Mb.
|
|
5.3. Спектральный и статистический анализ 209 Re{ ( )} Re{ ( )}; Im{ ( )} Im{ ( )} X m X m X m X m − = − = − , (5.13) т. е. действительная часть спектра является четной функцией частоты, а мни- мая часть спектра - нечетной функцией частоты. Разложения (5.12) и (5.8) позволяют рассматривать совокупность комплекс- ных амплитуд (5.9) как изображение периодического процесса в частотной облас- ти. Желание распространить такой подход на произвольные, в том числе - непе- риодические процессы привело к введению понятия Фурье-изображения в соот- ветствии со следующим выражением: X f x t e dt j f t ( ) ( ) ( ) = ⋅ − −∞ +∞ ∫ 2 π (5.14) Этот интеграл, несмотря на его внешнее сходство с выражением (5.9) для комплексных коэффициентов ряда Фурье, довольно существенно отличается от них. Во-первых, в то время, как физическая размерность комплексной амплиту- ды совпадает с размерностью самой физической величины x t ( ) , то размерность Фурье-изображения равна размерности x t ( ) , умноженной на размерность време- ни. Во-вторых, интеграл (5.14) существует (является сходящимся к конечной величине) только для так называемых "двухсторонне затухающих" процессов (т.е. таких, которые стремятся к нулю как при t → +∞ , так при t → −∞ ). Иначе гово- ря, его нельзя применять к так называемым "стационарным" колебаниям. Обратное преобразование Фурье-изображения в исходный процесс x t ( ) в этом случае определяется интегралом x t X f e df j f t ( ) ( ) ( ) = ⋅ −∞ +∞ ∫ 2 π , (5.15) который представляет собой некоторый аналог комплексного ряда Фурье (5.1). Указанное серьезное противоречие несколько сглаживается при численных расчетах, так как в этом случае можно иметь дело только с процессами ограни- ченной длительности, причем сам процесс в заданном диапазоне времени должен быть задан своими значениями в ограниченном числе точек. В этом случае интегрирование заменяется суммированием, и вместо вычис- ления интеграла (5.14) ограничиваются вычислением суммы X k f t x m t e j k m f m n [( ) ] [( ) ] ( ) ( ) − ⋅ = ⋅ − ⋅ ⋅ − ⋅ ⋅ − ⋅ − ⋅ ⋅ = ∑ 1 1 2 1 1 1 ∆ ∆ ∆ ∆ ∆ π t . (5.16) Тут, по сравнению с интегралом (5.14) осуществлены такие замены непрерывный интеграл приближенно заменен ограниченной суммой площадей прямоугольников, одна из сторон которых равна дискрету времени ∆t , с которым представлены значения процесса, а вторая - мгно- венному значению процесса в соответствующий момент времени; непрерывное время t заменено дискретными его значениями ( ) m t − ⋅ 1 ∆ , где - номер точки от начала процесса; m 5.3. Спектральный и статистический анализ 210 непрерывные значения частоты f заменены дискретными ее значениями ( ) k f − ⋅ 1 ∆ , где k - номер значения частоты, а дискрет частоты равен ∆f T = 1 , где T - промежуток времени, на котором задан процесс; дифференциал dt заменен ограниченным приращением времени ∆t . Если обозначить дискрет времени ∆t через Ts, ввести обозначения x m x m t ( ) [( ) = ] − ⋅ 1 ∆ ; X k X k f ( ) [( ) ] = − ⋅ 1 ∆ . а также учесть то, что число точек, в которых задан процесс, равно n n T t T Ts f t = = = ⋅ ∆ ∆ 1 ∆ 1 n , (5.17) то соотношение (5.15) можно представить в более удобной форме: X k Ts x m e j n k m m n ( ) ( ) ( / ) ( ) ( ) = ⋅ ⋅ − ⋅ ⋅ − ⋅ − = ∑ 2 1 1 π . (5.18) Как было отмечено в разделе 1.4.5 (формулы (2) и (3)), процедуры MatLAB fft и ifft осуществляют вычисления в соответствии с формулами: y k x m e j m k m n ( ) ( ) ( ) ( )/ = ⋅ − ⋅ ⋅ − ⋅ − = ∑ 2 1 1 1 π ; (5.19) x m n y k e j m k k n ( ) ( ) ( ) ( )/ = ⋅ ⋅ ⋅ − ⋅ − = ∑ 1 2 1 1 1 π n (5.20) соответственно. Сравнивая (5.18) с (5.19), можно сделать вывод, что процедура fft находит дискретное Фурье-изображение заданного дискретного во времени про- цесса x t ( ) , поделенное на дискрет времени: y k X k Ts ( ) ( ) = . (5.21) Осуществляя аналогичную операцию дискретизации соотношения (5.9) для комплексной амплитуды , получим X k * ( ) X k Ts T x m e j n k m m n * ( / ) ( ( ) ( ) = ⋅ ⋅ − ⋅ ⋅ − ⋅ − = ∑ 2 1 1 1 π ) ( ) = = ⋅ ⋅ = − ⋅ ⋅ − ⋅ − = ∑ 1 2 1 1 1 n x m e y k n j n k m m n ( ) ( ) ( / ) ( ) ( ) π (5.22) Из этого следует, что комплексный спектр разложения стационарного про- цесса равен поделенному на число измерений результату применения процедуры fft к заданному вектору измеренного процесса. Если же принять во внимание, что для большинства стационарных колеба- тельных процессов именно частотный, амплитудный и фазовый спектры не зави- сят от длительности T конкретной реализации и выбранного дискрета времени , то надо также сделать вывод, что Ts для спектрального анализа стационарных процессов наиболее целесообразно применять процедуру fft , результат которой делить затем на число точек измерений. 5.3. Спектральный и статистический анализ 211 Перейдем к определению Спектральной Плотности Мощности (СПМ), или, сокращенно, просто Спектральной Плотности (СП). Это понятие в теории опреде- ляется как Фурье-изображение так называемой корреляционной функции R 12 ( ) τ и применяется, в основном, для двух одновременно протекающих стационарных процессов и . Взаимная корреляционная функция (ВКФ) двух таких процессов определяется соотношением: x t 1 ( ) x t 2 ( ) R T x t x t dt T T T 12 1 2 2 2 1 ( ) lim ( ) ( ) τ = ⋅ →∞ − + ∫ τ + ⋅ , (5.23) т.е. ВКФ является средним во времени значением произведения первой функции на сдвинутую относительно нее на время задержки τ вторую функцию. Итак, Взаимная Спектральная Плотность (ВСП) двух стационарных процес- сов может быть определена так: S f R e d j f 12 12 2 ( ) ( ) ( ) = ⋅ − −∞ +∞ ∫ τ π τ τ (5.24) При числовых расчетах, когда оба процесса и заданы на опре- деленном ограниченном промежутке x t 1 ( ) x t 2 ( ) T времени своими значениями в некоторых точках, разделенных дискретом времени , формулу (5.23) можно трансфор- мировать в такую: n Ts R l n l x m x m l m n l 12 1 2 1 1 1 ( ) ( ) ( ) = − ⋅ + − = − ∑ , ( l n = 1 2 2 , , ... / ); (5.25) или в несколько более простое соотношение R l n x m x m l m n 12 1 2 1 2 2 1 ( ) ( ) ( ) / = ⋅ + = ∑ − , (l n = 1 2 2 , , ... / ); (5.26) а вместо (5.24) использовать S k Ts R l e j n k l l n 12 12 2 1 1 2 ( ) ( ) ( / ) ( ) ( ) / = ⋅ ⋅ − ⋅ ⋅ − ⋅ − = ∑ π 1 , ( k n = 1 2 2 , , ... / ). (5.27) Если теперь подставить выражение (5.26) в (5.27) и изменить в нем порядок суммирования, то можно прийти к такому соотношению между ВСП и результа- тами преобразований процедурой fft заданных измеренных значений процессов : S k Ts n y k y k 12 2 1 2 ( ) ( ) ( ) = ⋅ ⎧⎨ ⎩ ⎫ ⎬ ⎭ ⋅ ; ( k n = 1 2 2 , , ... / ), (5.28) где черта сверху означает комплексное сопряжение соответствующей величины. С учетом (5.21) и (5.22) выражение (5.28) можно представить также в виде: S k X k X k 12 2 1 ( ) ( ) ( ) * = ⋅ . (5.29) Из этого следует, что взаимная спектральная плотность двух процессов при любом значении частоты равна произведению значения комплексного спек- тра второго процесса на комплексно-сопряженное значение Фурье-изображения первого процесса на той же частоте. 5.3. Спектральный и статистический анализ 212 Формулы (5.21), (5.22) и (5.28) являются основой для вычислений в систе- ме MatLAB соответственно Фурье-изображения процесса, его комплексного спек- тра и взаимной спектральной плотности двух процессов. 5.3.2. Примеры спектрального анализа Чтобы применить процедуру fft как преобразование процесса, представлен- ного во временной области, в его представление в частотной области, следует, как было отмечено в разделе 1.4.5, сделать следующее: - по заданному значению дискрета времени Ts рассчитать величину F max диапазона частот (в герцах) по формуле: Fmax = 1/Ts ; (5.30) - по заданной длительности заданного процесса Т рассчитать дискрет частоты df по формуле: df =1/T ; (5.31) - по вычисленным данным сформировать вектор значений частот, в кото- рых будет вычислено Фурье-изображение. Последнее проще (но не наиболее правильно) сделать таким образом: f1=0 : df : Fmax. (5.32) В результате применения процедуры fft будет получено представление процесса в частотной области. Обратная процедура ifft , если ее применить к ре- зультатам первого преобразования, дает возможность восстановить исходный процесс во временной области. Однако процедура fft не дает непосредственно Фурье-изображения процес- са. Чтобы получить Фурье-изображение, надо сделать следующее (см. раздел 1.4.5): к результатам действия процедуры fft применить процедуру fftshift , которая переставляет местами первую и вторую половины полученного вектора; перестроить вектор частот по алгоритму f = -Fmax/2 : df : Fmax/2 . (5.33) Приведем примеры. Фурье-изображение прямоугольного импульса Сформируем процесс, состоящий из одиночного прямоугольного импульса. Зададим дискрет времени Ts=0.01с, длительность процесса Т=100с, амплитуду импульса А=0.75 и его ширину w=0.5с: Ts=0.01; T=100; A=0.75; w=0.5; t=0 : Ts : T; y = A*rectpuls(t, w); plot(t(1:100),y(1:100)), grid, set(gca,'FontName','Arial Cyr','FontSize',16), title('Процесс из одиночного прямоугольного импульса '); xlabel('Время (с)'); ylabel('Y(t)') Результат показан на рис. 5.27. 5.3. Спектральный и статистический анализ 213 Рис. 5.27 Рис. 5.28 Применим к вектору y процедуру fft и построим график зависимости мо- дуля результата от частоты. При этом графики в частотной области удобнее вы- водить при помощи процедуры stem (см. рис. 5.28): x=fft(y); df=1/T; Fmax=1/Ts; f=0 : df : Fmax; a=abs(x); stem(f,a), grid, set(gca,'FontName','Arial Cyr','FontSize',14), title('Модуль FFT-преобразования прямоугольного импульса '); xlabel('Частота (Гц)'); ylabel('Модуль') Теперь построим график модуля Фурье-изображения процесса: xp=fftshift(x); f1=-Fmax/2 : df : Fmax/2; a=abs(xp); stem(f1,a), grid, set(gca,'FontName','Arial Cyr','FontSize',14), title('Модуль Фурье-изображения прямоугольного импульса '); xlabel('Частота (Гц)'); ylabel('Модуль') Получим результат, приведенный на рис. 5.29. 5.3. Спектральный и статистический анализ 214 Рис. 5.29 Рис. 5.30 В заключение построим графики действительной и мнимой частей Фурье- изображения прямоугольного импульса: dch=real(xp); mch=imag(xp); plot(f1,dch,f1,mch), grid, set(gca,'FontName','Arial Cyr','FontSize',16), title('Фурье-изображение прямоугольного импульса '); ylabel('Действит. и Мнимая части'), xlabel('Частота (Гц)'); Они представлены на рис. 5.30. Фурье-изображение полигармонического процесса Рассмотрим пример трехчастотных гармонических колебаний - с частотою 1/ π , 1 та 3 Гц и амплитудами соответственно 0.6, 0.3 та 0.7: y t t t t ( ) . cos( ) . sin( ) . cos( / ) = ⋅ + ⋅ ⋅ + ⋅ ⋅ + 0 6 2 0 3 2 0 7 6 4 π π π 5.3. Спектральный и статистический анализ 215 Найдем Фурье-изображение этого процесса и выведем графики самого про- цесса, модуля его Фурье-изображения, а также действительную и мнимую части: Ts = 0.01; T = 100; t = 0 : Ts : T; Y = 0.6*cos(2*t)+0.3*sin(2*pi*t)+0.7*cos(6*pi*t+pi/4); plot(t,Y), grid, set(gca,'FontName','Arial Cyr','FontSize',16), title(' Трехчастотный полигармонический процесс '); xlabel('Время (с)'); ylabel('Y(t)') График процесса показан на рис. 5.31. Рис. 5.31 Находим модуль Фурье-изображения этого процесса: df = 1/T; Fmax = 1/Ts; dovg=length(t); f = - Fmax/2 : df : Fmax/2; X = fft(Y); Xp = fftshift(X); A = abs(Xp); s1 = dovg/2 - 400; s2 = dovg/2 + 400; stem(f(s1:s2),A(s1:s2)), grid, set(gca,'FontName','Arial Cyr','FontSize',14), title('Модуль Фурье-изображения полигармонического процесса'); xlabel('Частота (Гц)'); ylabel('Модуль') Результат представлен на рис. 5.32. 5.3. Спектральный и статистический анализ 216 Рис. 5.32 Если изменить дискрет времени на Ts=0.02, получим результат, изобра- женный на рис. 5.33. Рис. 5.33 Как видно, результат Фурье-преобразования в значительной степени зави- сит от величины дискрета времени и мало что говорит об амплитудах гармониче- ских составляющих. Это обусловлено различием между определениями Фурье- изображения и комплексного спектра. Поэтому для незатухающих (установив- шихся, стационарных) колебаний любого вида намного удобнее находить не Фу- рье-изображение, а его величину, деленную на число точек в реализации. В пре- дыдущей части программы это эквивалентно замене оператора X=fft(Y) на X = fft(Y)/dovg , где dovg - длина вектора t. 5.3. Спектральный и статистический анализ 217 Рис. 5.34 Рис. 5.35 В результате получается комплексный спектр (рис. 5.34), полностью соот- ветствующий коэффициентам комплексного ряда Фурье. Выделим действительную и мнимую части комплексного спектра: dch = real(Xp); mch = imag(Xp); s1 = dovg/2 - 400; s2 = dovg/2 + 400; subplot(2,1,1), plot(f(s1:s2),dch(s1:s2)), grid, set(gca,'FontName','Arial Cyr','FontSize',10), title('Комплексный спектр полигармонических колебаний'); ylabel('Действит. часть'); subplot(2,1,2) , plot(f(s1:s2),mch(s1:s2)), grid, set(gca,'FontName','Arial Cyr','FontSize',10), 5.3. Спектральный и статистический анализ 218 xlabel('Частота (Гц)'); ylabel('Мнимая часть') По полученным графикам (рис. 5.35) можно судить не только о частотах и амплитудах, а и о начальных фазах отдельных гармонических составляющих. Фурье-изображение случайного процесса В заключение рассмотрим Фурье-преобразование случайного стационарно- го процесса, сформированного ранее (см. рис. 5.25). Сначала сформируем процесс в виде белого гауссового шума (рис. 5.36) с шагом во времени 0.01 и длительно- стью в 100 с : Ts=0.01; T = 100; t=0 : Ts : T; % Задание параметров процесса x1=randn(1,length(t)); % Формирование белого шума % Построение графика белого шума plot(t,x1),grid, set(gca,'FontName','Arial Cyr','FontSize',16) title('Белый Гауссовый шум (СКО= 1; Ts= 0.01)'); |