logo
Цифровая обработка сигналов

5.3.2. Примеры спектрального анализа

Чтобы применить процедуру fft как преобразование процесса, представленного во временной области, в его представлении в частотной области, необходимо, как было отмечено в разделе 1.4.5., сделать следующее:

(5.30)

(5.31)

Последнее проще (но не наиболее правильно) сделать таким образом:

(5.32)

В результате применения процедуры fft будет получено представление процесса в частотной области. Обратная процедура ifft, если ее применить к результатам первого преобразования, дает возможность восстановить исходный процесс во временной области.

Однако процедура fft не дает непосредственно Фурье-изображения процесса. Чтобы получить Фурье-изображение, необходимо выполнить следующие действия (разд. 1.4.5):

(5.33)

Сформируем процесс, состоящий из одиночного прямоугольного импульса. Зададим дискрет времени Ts=0.01 c, длительность процесса Т=100 с, амплитуду импульса А=0.75 и его ширину w=0/5 c:

Ts=0.01; T=100; A=0.75; w=0.5;

t=o : 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. /

Применим к вектору 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 (`Модуль`)

Теперь построим график модуля Фурье-изображения процесса?

x = 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.

В заключении построим графики действительной и мнимой частей Фурье-изображения прямоугольного импульса:

dch=real (xp); mch=imag(xp);

plot(f1, dch, mch), grid, set (gca, `FontName`, …

`Arial Cyr`, `FontSize`, 16),

title (`Фурье-изображениt прямоугольного импульса`);

ylabel (`Действит. и Мнимая части`)

xlabel (`Частота (Гц)`)

Полученные графики представлены на рис. 5.30.

Фурье-изображение полигармонического процесса

Рассмотрим пример трехчастотных гармонических колебаний – с частотой , 1 и 3 Гц и амплитудами соответственно 0.6, 0.3 и 0.7:

Найдем Фурье-изображение этого процесса и выведем графики самого процесса, модуля его Фурье-изображения, а также действительную и мнимую части:

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.

Находим модуль Фурье-изображения этого процесса:

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.

Если изменить дискрет времени на Ts=0.02, получим результат, изображенный на рис. 5.33.

Как видно, результат Фурье-преобразования в значительной степени зависит от величины дискрета времени и мало что говорит об амплитудах гармонических составляющих. Это обусловлено различием между определениями Фурье-изображения и комплексного спектра. Поэтому для незатухающих (установившихся, стационарных) колебаний любого вида намного удобнее находить не Фурье-изображение, а его величину, деленную на число точек в реализации. В предыдущей части программы это эквивалентно замене оператора X=fft(Y) на X=fft(Y)/dovg, где dovg – длина вектора t.

В результате получается комплексный спектр (рис. 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),

xlabel (`Частота (Гц)`)

ylabel (`Мнимая часть`)

По полученным графикам (рис. 5.35.) можно судить не только о частотах и амплитудах, но и о начальных фазах отдельных гармонических составляющих.

Фурье-изображение случайного процесса

В заключении рассмотрим Фурье-преобразование случайного стационарного процесса, сформированного ранее (рис. 5.25). Сначала сформируем процесс в виде белого гауссового шума (рис. 5.36) с шагом во времени 0.01 и длительностью 100 с:

Ts=0.01; T=100; % Задание параметров процесса

t=0 : Ts : T;

x1=randn(l, length(t); % Формирование белого шума

% Построение графика белого шума

plot(t, x), grid, set (gca, `FontName`, `Arial Cyr`, `FontSize`, 16),

title (`Белый Гауссовый шум (СКО=1ж; Ts=0.01)`)

xlabel (`Время (с)`)

ylabel (`X1(t)`);

Теперь создадим формирующий фильтр, «пропустим» через него белый шум и результат выведем на рис. 5.37:

% Расчет параметров формирующего фильтра

omO=2*pi; dz=0.05; A1=oms=omO*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)

titlt (`Процесс на выходе фильтра (TO=1; dz=0.05, Ts=0.01)`);

xlabel (`Время (с)`);

ylabel (`Y1(t)`)

Вычислим Фурье-изображение (ФИ) для процесса-шума с учетом замечания, сделанного для установившихся процессов, и построим на рис. 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 (`Модуль ФИ гауссового белого шума`);

subplot(2, 1, 2);

stem(f, S1), grid,

set (gca, `FontName`, `Arial Cyr`, `FontSize`, 10)

title (`Спектральная плотность мощности`);

xlabel(`Частота (Гц)`)

Из рис. 5.38 видно, что спектральная плотность практически одинакова по величине во всем диапазоне частот, чем и обусловлено название процесса – белый шум.

Аналогичную процедуру выполним и для «профильтрованного» процесса (рис. 5.39):

% Вывод графиков профильтрованного процесса

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)

title (`Спектральная плотность мощности`);

xlabel(`Частота (Гц)`)

Проводя эти вычисления еще раз с новой длительностью процесса Т=20 с (рис. 5.40), можно наглядно убедиться что величины ФИ и СПМ при этом практически не изменяются.