Построение БПФ на октаве

Я знаю, что FFT изменяет функцию во временной области на одну, показанную в частотной области.

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

Кроме того, я могу только получить амплитуды, чтобы они соответствовали значениям в исходном сигнале, разделив ось y на некоторое целое число. Почему это?

Здесь мой код

t=0:0.001:2

x=2*sin(20*pi*t) + sin(100*pi*t)
subplot(2,1,1)
plot(1000*t,x)
grid
xlabel("Time in milliseconds")
ylabel("Signal amplitude")

subplot(2,1,2)
y=fft(x)
plot(1000*t,abs(y))
xlabel("Frequency")
ylabel("Signal amplitude")

и графики.

enter image description here

Пожалуйста, помогите = (

Ответы

Ответ 1

Отношение частоты (масштабирование по оси x)

Частота каждого значения, полученного БПФ, линейно связана с индексом выходного значения через:

f(i) = (i-1)*sampling_frequency/N

Где N - число точек FFT (т.е. N=length(y)). В вашем случае N=2001.

Можно вычесть частоту дискретизации из вашего определения t как 1/T, где T - временной интервал выборки (T = 0,001 в вашем случае). Таким образом, частота дискретизации составляет 1000 Гц.

Заметим, что поскольку значение t(i) также линейно связано с индексом i, через

t(i) = (i-1)*0.001

можно (хотя не обязательно советовать, так как это просто заслонит ваш код), чтобы определить f = 1000*t*sampling_frequency/N. Обратите внимание, что вам не хватает члена sampling_frequency/N, который соответственно приводил к отображению тонов с неправильной частотой (из определения x должны быть пики с частотой 10 Гц и 50 Гц и соответствующие алиасы на частотах 990 Гц и 950 Гц).

Отношение амплитуды (масштабирование по оси Y)

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

Упрощение проблемы до одного тона:

x = A*sin(2*pi*f*t)

Приблизительная амплитуда соответствующего пика может быть получена с помощью Теорема Парсеваля:

enter image description here

Во временной области (левая часть уравнения) выражение приблизительно равно 0.5*N*(A^2).

В частотной области (правая часть уравнения), делая следующие предположения:

  • эффекты спектральной утечки являются незначительными.
  • спектральное содержание тона содержится всего в 2 булавках (на частоте f и соответствующей частоте сглаживания sampling_frequency-f) учитывается суммирование (все остальные бины равны ~ 0). Обратите внимание, что это обычно выполняется только в том случае, если частота тона является точным (или почти точным) кратным sampling_frequency/N.

выражение в правой части приблизительно равно 2*(1/N)*abs(X(k))^2 для некоторого значения k, соответствующего пику на частоте f.

Полагая два вместе, получаем abs(X(k)) ~ 0.5*A*N. Другими словами, выходная амплитуда показывает коэффициент масштабирования 0.5*N (или приблизительно 1000 в вашем случае) относительно амплитуды во временной области, как вы это заметили.

Идея по-прежнему применяется с более чем одним тоном (хотя пренебрежимо малое предположение о спектральной утечке в конечном итоге ломается).

Ответ 2

В других ответах было высказано предположение, что в этом примере имеются частотные характеристики на частотах 950 Гц и 990 Гц. Это недоразумение о том, как код FFT использует индексы. Эти "высокочастотные" пики на самом деле -50Hz и -10Hz.

Частотная область простирается от -N/2 * sampling_frequency/N до + N/2 * sampling_frequency/N. Но по историческим причинам принято, что первые N/2 фрагментов информации являются положительными частотами, средняя точка - нулевой частотой, а последние N/2 фрагментов информации - отрицательными частотами в обратном порядке. Для спектра мощности нет необходимости показывать больше, чем первые 1 + N/2 фрагмента информации.

Это соглашение чрезвычайно запутанно, так как мне пришлось его озадачить из Press et al. Численные рецепты и ручное кодирование Быстрого Преобразования Хартли, много лет назад, когда я впервые использовал БПФ, еще до выпуска бета-теста Matlab 1.0, который Клив Молер раздал некоторым счастливчикам :-)

Ответ 3

Вы можете использовать fftshift для центровки с нулевой частотой, но это не так важно.

Согласно теореме Найквиста/Котельникова, самая высокая частота, которую вы получаете, составляет половину частоты дискретизации: fmax = fsample/2. Поскольку у вас есть частота дискретизации 1 кГц, вы получаете максимум 500 Гц. Шаг частоты - это инверсия времени приема: df = 1/2 = 0,5 Гц в вашем случае.

Итак, формула, указанная в другом ответе, неверна, вы должны ее вдвое.

Теперь у вас 1000 посередине графика. Вы будете использовать только первую половину. Для каждой точки вы получаете половину номера на оси x, чтобы получить частоту.

PS.I не мог комментировать, потому что у меня недостаточно репутации.