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

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

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

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

$ω = e^{- 2 π i / n}$ один из $n$ комплексные корни из единицы, где $i$ мнимая единица. Для $x$ и $y$ , индексы $j$ и $k$ расположитесь от $0$ к $n - 1$ .

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

t = 0:1/50:10-1/50;                     
x = sin(2*pi*15*t) + sin(2*pi*20*t);
plot(t,x)

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

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

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

plot(f,abs(y))
title('Magnitude')

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

n = length(x);                         
fshift = (-n/2:n/2-1)*(50/n);
yshift = fftshift(y);
plot(fshift,abs(yshift))

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

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

xnoise = x + 2.5*gallery('normaldata',size(t),4);

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

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

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

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

Считайте аудиоданные собранными из подводных микрофонов недалеко от берегов Калифорнии. Эти данные могут быть найдены в библиотеке, сохраняемой Программой исследований Биоакустики Корнелльского университета. Загрузите и отформатируйте подмножество данных в bluewhale.au, который содержит Тихоокеанскую вокализацию голубого кита. Можно использовать команду 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)])

Задайте новую длину сигнала, которая является следующей степенью 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')
ylabel('Power')

Документация