Преобразование Фурье

Открыть пример в Engee

Преобразование Фурье – это формула, при помощи которой выполняется перевод одного вектора чисел в другой. Первый вектор состоит из набора значений сигнала, взятых в отдельные моменты времени. Второй вектор состоит из дискретного набора значений, характеризующих амплитуду и частоту гармонических сигналов, из которых можно составить наблюдаемый сигнал путем его реконструкции во времени или в пространстве.

Рассмотрим вектор $x$, включающий $n$ измерений, взятых через одинаковые интервалы времени. Преобразование Фурье такого вектора описано следующим образом:

$$y_{k+1} = \sum_{j=0}^{n-1}\omega^{jk}x_{j+1}$$

Угловая частота $\omega = e^{-2\pi i/n}$ является комплексным корнем из единицы, $i$ – мнимая единица. Для обоих векторов $x$ и $y$, индексы $j$ и $k$ принимают значения от 0 до $n-1$.

Функция fft из библиотеки FFTW выполняет разложение вектора в ряд Фурье при помощи алгоритма быстрого преобразования Фурье (БПФ). Рассмотрим синусоидальный сигнал $x$, каждому взятый через равные интервалы времени. Моменты осуществления измерений образуют вектор $t$. Сигнал состоит из двух компонентов с частотами 15 и 20 Гц. Изучим сигнал, измерение которого производилось каждую 1/50 долю секунды в течение 10 секунд.

Ts = 1/50;
t = 0:Ts:10-Ts;                     
x = sin.(2pi * 15t) .+ sin.( 2pi * 20t );
plot( t, x, xlabel="Время (с)", ylabel="Амплитуда", leg=false )

Разложим сигнал в ряд Фурье и создадим вектор $f$, дискретные компоненты которого образуют выборку сигнала в частотной области.

using FFTW

y = fft(x);
fs = 1/Ts;
f = (0:length(y)-1)*fs/length(y);

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

plot( f, abs.(y), title="Амплитудный спектр", xlabel="Частота (Hz)", ylabel="Амплитуда" )

На графике видны четыре всплеска, хотя исходный сигнал имел две составляющих с частотами 15 и 20 Гц. Дело в том, что правая половина графика зеркально отражает левую половину. Дискретное преобразование Фурье сигнала во временной области тоже имеет свою периодичность – первая половина вектора соответствует области "положительных значений частоты", а вторая половина – области "отрицательных значений частоты" составляющих. Всплески в области 30 и 35 Гц на самом деле соответствуют значениям -20 и -15 Гц.

На практике часто достаточно использовать только одну половину ряда Фурье, не выполняя лишних расчетов. Для этого можно пользоватьcя функцией rfft. Но для наглядности выведем весь ряд Фурье, сместив правую половину в область отрицательных частот при помощи функции fftshift.

n = length(x);
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot( fshift, abs.(yshift), title="Амплитудный спектр", xlabel="Частота (Uw)", ylabel="Амплитуда" )

Зашумленный сигнал

В реальном сценарии, измеренный сигнал часто включает как полезную составляющую, так и случайный шум, который не позволяет однозначно различить гармонические составляющие. Использование преобразования Фурье может помочь отсечь случайную составляющую и оставить только систематический, колебательный процесс. Рассмотрим новый вектор измерений xnoise, результат сложения сигнала x и случайного гауссовского шума.

using Random
Random.seed!(0)
xnoise = x .+ 2.5 .* randn( size(t) );

На практике, зависимость мощности сигнала от частоты используется чаще, чем зависимость амплитуды от частоты. Мощность равняется квадрату амплитуды компонента ряда Фурье, разделенная на количество точек в ряде. Рассчитаем и построим график спектра мощности сигнала, с нулевой частотой по центру оси $x$. Несмотря на наличие шума, мы ясно видим оба всплеска на графике мощности сигнала.

ynoise = fft( xnoise );
ynoiseshift = fftshift( ynoise );    
power = abs.( ynoiseshift ).^2 ./ n; 
plot( fshift, power, title="Спектр мощности", xlabel="Частота (Гц)", ylabel="Мощность сигнала" )

Скорость преобразования

Чтобы рассчитать преобразование Фурье всего ряда из $n$ элементов вектора $y$ при помощи исходной формулы, нужно выполнить порядка $n^2$ вычислений с числами с плавающей запятой. Алгоритм быстрого преобразования Фурье требует выполнить порядка $n log n$ действий, что делает этот метод очень привлекательным из-за экономии времени, особенно если обрабатываемый сигнал состоит из нескольких миллионов точек. Более того, конкретные реализации алгоритма БПФ могут вычисляться еще быстрее в ряде конкретных случаем. Например, в случае, когда $n$ равно 2 в некоторой степени.

Рассмотрим аудиофайл bluewhale.wav – запись акустических коммуникаций синих китов, записанных на подводный микрофон.

using WAV
x, fs = wavread( "bluewhale.wav" );

Синие киты общаются в инфразвуке, полоса частот акустических колебаний почти не пересекается с частотами, к которым чувствителен человеческий слух. Поэтому вектор времени в представленных данных сжат в 10 раз, чтобы высота звукового сигнала была доступна нашему восприятию.

Следующая ячейка позволяет воспроизвести звуковой сигнал, прочитанный из файла.

using Base64

buf = IOBuffer();
wavwrite(x, buf; Fs=fs);
data = base64encode(unsafe_string(pointer(buf.data), buf.size));
markup = """<audio controls="controls" {autoplay}>
              <source src="data:audio/wav;base64,$data" type="audio/wav" />
              Your browser does not support the audio element.
              </audio>"""
display( "text/html", markup );

Визуализируем этот сигнал:

whaleMoan = x[24500:31000];
t = 10 * (0:1/fs:(length(whaleMoan))/fs);
plot( t, whaleMoan, xlabel="Время (с)", ylabel="Амплитуда", label=false, xlimits=(0,t[end]) )

Теперь увеличим длину сигнала до ближайшего числа $2^N$ в сторону возрастания. Мы получим более длинный отрезок сигнала. И рассчитаем ряд Фурье для этого нового отрезка сигнала. Если бы длина сигнала не была взята из ряда $2^N$ (при $N = 1,2\dots$), то функция fft, скорее всего, все равно дополнила бы сигнал нулями до ближайшей длины $2^N$. Это делается для того, чтобы существенно ускорить алгоритм БПФ, особенно если значение количества точек не имеет малых делителей.

m = length( whaleMoan ); 
n = nextpow(2, m);

whaleMoan_zeros = zeros( n )
whaleMoan_zeros[1:length(whaleMoan)] = whaleMoan

y = fft( whaleMoan_zeros );

Построим график спектра мощности этого сигнала. Согласно графику, акустический сигнал включает основную гармонику – звук на частоте примерно 17 Гц, и несколько других заметных гармонических составляющих.

f = (0:n-1) .* (fs/n)/10; # вектор значений частоты
power = abs.(y).^2/n;     # спектр мощности
plot( f[1:Int(floor(n/2))],power[1:Int(floor(n/2))],
    xlabel = "Частота (Гц)", ylabel = "Мощность", leg=false )

Определение фазы синусоид

При помощи преобразования Фурье можно также построить фазовый спектр изучаемого сигнала. Ради примера мы создадим сигнал, в котором суммируется два гармонических сигнала с частотой 15 и 20 Гц. Первый сигнал будет косинусом со значением фазы $-\pi/4$, а второй сигнал будет синусоидой без фазового смещения. Частота дискретизации сигнала будет 100 Гц, длина сигнала 1 с.

fs = 100;
t = 0:1/fs:1-1/fs;
x = @. cos(2*pi*15*t - pi/4) - sin(2*pi*40*t);

Разложим сигнал в ряд Фурье и построим график амплитуды в зависимости от частоты.

y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*fs;

plot( f, abs.(z), st=:stem, marker=(5, :white, :o),
    markerstrokecolor=1,
    xlabel="Частота (Гц)", ylabel="|y|", grid=true )

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

tol = 1e-6;
z[abs.(z) .< tol] .= 0;

theta = angle.(z);

plot( f, theta/pi, st=:stem, marker=(5, :white, :o),
    markerstrokecolor=1,
    xlabel="Частота (Гц)", ylabel = "Фаза / π", grid=true )

Заключение

Мы построили несколько разных спектров (амплитуды, мощности и фазы) на основе разложения сигнала в ряд Фурье. В этом нам помогла функция fft, использующая алгоритм БПФ, и несколько других функций (например fftshift). В итоге мы смогли наглядно показать все интересующие нас свойства входного сигнала.