Как правильно масштабировать частотную ось в быстром преобразовании Фурье?

Question

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

import numpy as np
from matplotlib import pyplot as plt

N = 1024
limit = 10
x = np.linspace(-limit, limit, N)
dx = x[1] - x[0]
y = np.sin(2 * np.pi * 5 * x) + np.sin(2 * np.pi * x)
Y = np.abs(np.fft.fft(y) ** 2)
z = fft.fftshift(np.fft.fftfreq(N, dx))
plt.plot(z[int(N/2):], Y[int(N/2):])
plt.show()

Из заданной функции ясно, что должно быть два пика на частотах 1 и 5. Однако, когда я запускаю этот код, я получаю следующий график.

Очевидно, что шипы не там, где они должны быть. Кроме того, я заметил, что масштабирование частоты чувствительно к количеству точек N, а также к пределам интервала, которые я устанавливаю limit. Например, установка N = 2048 дает следующий график.

Как видите, расположение шипов изменилось. Теперь сохранение N = 1024 и установка limit = 100 также изменяет результат.

Как я могу сделать так, чтобы ось частоты все время корректно масштабировалась?

Answer 1

fftfreq возвращает частотный диапазон в следующем порядке: положительные частоты от самой низкой до самой высокой, затем отрицательные частоты в обратном порядке абсолютного значения. (Как правило, вы хотите построить только половину, как в своем коде.) Обратите внимание, что на самом деле функция должна очень мало знать о данных: только количество выборок и их интервал во временной области.

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

Однако он принимает комплексные числа в качестве входных данных. На практике это редко. Входные данные обычно представляют собой образцы действительных чисел, как в примере выше. В этом случае преобразование Фурье обладает особым свойством: оно симметрично в частотной области, т. Е. Имеет одинаковое значение для f и −f. По этой причине часто не имеет смысла отображать обе половины спектра, поскольку они содержат одинаковую информацию.

Есть одна выделяющаяся частота: f = 0. Это мера среднего значения сигнала, его смещение от нуля. В спектре, возвращаемом fft, и частотном диапазоне от fftfreq, он находится в самом первом индексе массива. Если строит графики для обеих половин, может иметь смысл сместить частотный спектр так, чтобы отрицательная половина была слева от нулевого компонента, положительная половина - справа, то есть все значения находятся в порядке возрастания заказать и готов к участку.

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

На самом деле вам не нужно полагаться ни на одну из этих функций, работающих на частотах. Сгенерировать их диапазон легко на основе документации только fftfreq:

from numpy.fft import fft
from numpy import arange, linspace, sin, pi as π
from matplotlib import pyplot

def FFT(t, y):
    n = len(t)
    Δ = (max(t) - min(t)) / (n-1)
    k = int(n/2)
    f = arange(k) / (n*Δ)
    Y = abs(fft(y))[:k]
    return (f, Y)

t = linspace(-10, +10, num=1024)
y = sin(2*π * 5*t) + sin(2*π * t)
(f, Y) = FFT(t, y)
pyplot.plot(f, Y)
pyplot.show()

Обратите внимание, что Numpy также предлагает специальные функции, rfft rfftfreq, для общего случая использования реальных данных.