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

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

Преобразование Фурье задано для вектора x с n равномерно выбранные точки

yk+1=j=0n-1ωjkxj+1.

ω=e-2πi/n является одним из n комплексные корни единства, где i - мнимый модуль. Для x и y, индексы j и k область значений от 0 кому n-1.

The fft функция в MATLAB ® использует быстрый алгоритм преобразования Фурье, чтобы вычислить преобразование Фурье данных. Рассмотрим синусоидальный сигнал x это функция времени t с частотными составляющими 15 Гц и 20 Гц. Используйте временной вектор, дискретизированный с шагами 150 секунды в течение 10 секунд.

Ts = 1/50;
t = 0:Ts:10-Ts;                     
x = sin(2*pi*15*t) + sin(2*pi*20*t);
plot(t,x)
xlabel('Time (seconds)')
ylabel('Amplitude')

Figure contains an axes. The axes contains an object of type line.

Вычислите преобразование Фурье сигнала и создайте вектор f что соответствует дискретизации сигнала в разносе частот.

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

Когда вы строите график величины сигнала как функции от частоты, пики в величине соответствуют частотным составляющим сигнала 15 Гц и 20 Гц.

plot(f,abs(y))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Magnitude')

Figure contains an axes. The axes with title Magnitude contains an object of type line.

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

n = length(x);                         
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot(fshift,abs(yshift))
xlabel('Frequency (Hz)')
ylabel('Magnitude')

Figure contains an axes. The axes contains an object of type line.

Сигналы с шумом

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

rng('default')
xnoise = x + 2.5*randn(size(t));

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

ynoise = fft(xnoise);
ynoiseshift = fftshift(ynoise);    
power = abs(ynoiseshift).^2/n; 
plot(fshift,power)
title('Power')
xlabel('Frequency (Hz)')
ylabel('Power')

Figure contains an axes. The axes with title Power contains an object of type line.

Вычислительная эффективность

Использование формулы преобразования Фурье непосредственно для вычисления каждой из n элементы y требует от порядка n2 операции с плавающей точкой. Быстрое преобразование Фурье требует только порядка nlogn операции для вычисления. Эта вычислительная эффективность является большим преимуществом при обработке данных, которые имеют миллионы точек данных. Многие специализированные реализации алгоритма быстрого преобразования Фурье еще более эффективны, когда n является степенью 2.

Рассмотрим аудио данных, собранные с подводных микрофонов у побережья Калифорнии. Эти данные можно найти в библиотеке, поддерживаемой Программой исследований биоакустики Корнеллского университета. Загрузка и форматирование подмножества данных в bluewhale.au, который содержит тихоокеанскую вокализацию голубых китов. Поскольку вызовы синих китов являются низкочастотными звуками, они едва слышны для людей. Шкала времени в данных сжимается в 10 раз, чтобы поднять тангаж и сделать вызов более четко слышимым. Можно использовать команду sound(x,fs) прослушивание всего аудио файла.

whaleFile = 'bluewhale.au';
[x,fs] = audioread(whaleFile);
whaleMoan = x(2.45e4:3.10e4);
t = 10*(0:1/fs:(length(whaleMoan)-1)/fs);

plot(t,whaleMoan)
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 t(end)])

Figure contains an axes. The axes contains an object of type line.

Задайте новую длину сигнала, которая является следующей степенью 2, больше исходной длины. Затем используйте fft вычислить преобразование Фурье с помощью новой длины сигнала. fft автоматически заполняет данные нулями, чтобы увеличить размер выборки. Это заполнение может сделать расчет преобразования значительно быстрее, особенно для размеров выборки с большими простыми множителями.

m = length(whaleMoan); 
n = pow2(nextpow2(m));
y = fft(whaleMoan,n);        

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

f = (0:n-1)*(fs/n)/10; % frequency vector
power = abs(y).^2/n;   % power spectrum      
plot(f(1:floor(n/2)),power(1:floor(n/2)))
xlabel('Frequency (Hz)')
ylabel('Power')

Figure contains an axes. The axes contains an object of type line.

См. также

| | | | | |

Похожие темы