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

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

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

Угловая частота является комплексным корнем из единицы, – мнимая единица. Для обоих векторов и , индексы и принимают значения от 0 до .

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

Pkg.add(["WAV"])

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

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

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) );

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

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

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

Чтобы рассчитать преобразование Фурье всего ряда из элементов вектора при помощи исходной формулы, нужно выполнить порядка вычислений с числами с плавающей запятой. Алгоритм быстрого преобразования Фурье требует выполнить порядка действий, что делает этот метод очень привлекательным из-за экономии времени, особенно если обрабатываемый сигнал состоит из нескольких миллионов точек. Более того, конкретные реализации алгоритма БПФ могут вычисляться еще быстрее в ряде конкретных случаем. Например, в случае, когда равно 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]) )

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

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 Гц. Первый сигнал будет косинусом со значением фазы , а второй сигнал будет синусоидой без фазового смещения. Частота дискретизации сигнала будет 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). В итоге мы смогли наглядно показать все интересующие нас свойства входного сигнала.