Как рассчитать частоту данных с помощью БПФ? - PullRequest
Новая
       10

Как рассчитать частоту данных с помощью БПФ?

8 голосов
/ 19 ноября 2010

Я хочу знать частоту данных. Я немного подумал, что это можно сделать с помощью БПФ, но я не уверен, как это сделать. Как только я передал полные данные в FFT, он дает мне 2 пика, но как я могу получить частоту?

Заранее большое спасибо.

Ответы [ 6 ]

10 голосов
/ 20 ноября 2010

Предположим, x[n] = cos(2*pi*f0*n/fs), где f0 - частота вашей синусоиды в герцах, n=0:N-1, а fs - частота дискретизации x в выборках в секунду.

Пусть X = fft(x). И x, и X имеют длину N. Предположим, что X имеет два пика в n0 и N-n0.

Тогда синусоидальная частота равна f0 = fs*n0/N Гц.

Пример : fs = 8000 выборок в секунду, N = 16000 выборок. Следовательно, x длится две секунды.

Предположим, что X = fft(x) имеет пики на 2000 и 14000 (= 16000-2000). Следовательно, f0 = 8000 * 2000/16000 = 1000 Гц.

8 голосов
/ 09 июня 2011

Вот то, что вы, вероятно, ищете:

Когда вы говорите о вычислении частоты сигнала, вы, вероятно, не так заинтересованы в составляющих синусоидальных волнах. Это то, что дает вам БПФ. Например, если вы суммируете sin (2 * pi * 10x) + sin (2 * pi * 15x) + sin (2 * pi * 20x) + sin (2 * pi * 25x), вы, вероятно, захотите определить «частоту» "как 5 (взгляните на график этой функции). Однако БПФ этого сигнала будет обнаруживать величину 0 для частоты 5.

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

  1. Окно сигнала, используя гладкое окно (косинус подойдет. Окно должно быть как минимум в два раза больше, чем самый большой период, который вы хотите обнаружить. 3 раза больше даст лучшие результаты). (см. http://zone.ni.com/devzone/cda/tut/p/id/4844, если вы не уверены).

  2. Возьмите БПФ (однако, убедитесь, что размер БПФ в два раза больше окна, а вторая половина заполнена нулями. Если размер БПФ равен только размеру окна, вы будете эффективно беру круговую автокорреляцию, а это не то, что вам нужно. см. https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Circular_convolution_theorem_and_cross-correlation_theorem)

  3. Заменить все коэффициенты БПФ на их квадратные значения (действительные ^ 2 + imag ^ 2). Это эффективно принимает автокорреляцию.

  4. Возьми iFFT

  5. Найдите самый большой пик в iFFT. Это самая сильная периодичность сигнала. На самом деле вы можете быть немного более умным в том, какой пик вы выбираете, но для большинства целей этого должно быть достаточно. Чтобы найти частоту, достаточно взять f = 1 / T.

6 голосов
/ 01 апреля 2012

Если у вас есть сигнал с одной частотой (например: 

y = sin(2 pi f t)

С:

  • y сигнал времени
  • f центральной частоты
  • t time

Тогда вы получите два пика, один с частотой, соответствующей f, и один с частотой, соответствующей -f.

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

В коде: (я написал это на python, но он должен быть одинаково прост в c #): 

import numpy as np
from pylab import * 
x = np.random.rand(100) # create 100 random numbers of which we want the fourier transform
x = x - mean(x) # make sure the average is zero, so we don't get a huge DC offset.
dt = 0.1 #[s] 1/the sampling rate
fftx = np.fft.fft(x) # the frequency transformed part
# now discard anything  that we do not need..
fftx = fftx[range(int(len(fftx)/2))]
# now create the frequency axis: it runs from 0 to the sampling rate /2
freq_fftx = np.linspace(0,2/dt,len(fftx))
# and plot a power spectrum
plot(freq_fftx,abs(fftx)**2)
show()

Теперь частота находится на самом большом пике.

4 голосов
/ 20 ноября 2010

Если вы смотрите на величину, полученную из БПФ наиболее часто используемого типа, то сильная синусоидальная частотная составляющая реальных данных будет отображаться в двух местах, один раз в нижней половине, плюс его комплексное сопряженное зеркальное изображение в верхняя половина. Эти два пика представляют один и тот же спектральный пик и одну и ту же частоту (для строго реальных данных). Если номера бина результатов FFT начинаются с 0 (ноль), то наиболее вероятна частота синусоидальной составляющей, представленной бином в нижней половине результата FFT. 

Frequency_of_Peak = Data_Sample_Rate * Bin_number_of_Peak / Length_of_FFT ;

Убедитесь, что вы отработали правильные единицы в приведенном выше уравнении (чтобы получить единицы циклов в секунду, за две недели, за килопарсек и т. Д.)

Обратите внимание, что, если длина волны данных не равна точному целочисленному кратному длине БПФ, фактический пик будет между бинами, таким образом распределяя энергию между несколькими соседними бинами результата БПФ. Поэтому вам, возможно, придется интерполировать, чтобы лучше оценить пик частоты. Обычными методами интерполяции для нахождения более точной оценки частоты являются 3-точечная параболическая и свертка Синка (что почти совпадает с использованием более длинного БПФ с нулевым дополнением).

1 голос
/ 24 февраля 2014

Предполагая, что вы используете дискретное преобразование Фурье для просмотра частот, вы должны быть осторожны с тем, как интерпретировать нормализованные частоты обратно в физические (т. Е. Гц).

В соответствии с руководством FFTW о том, как рассчитать спектр мощности сигнала: 

#include <rfftw.h>
...
{
     fftw_real in[N], out[N], power_spectrum[N/2+1];
     rfftw_plan p;
     int k;
     ...
     p = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
     ...
     rfftw_one(p, in, out);
     power_spectrum[0] = out[0]*out[0];  /* DC component */
     for (k = 1; k < (N+1)/2; ++k)  /* (k < N/2 rounded up) */
          power_spectrum[k] = out[k]*out[k] + out[N-k]*out[N-k];
     if (N % 2 == 0) /* N is even */
          power_spectrum[N/2] = out[N/2]*out[N/2];  /* Nyquist freq. */
     ...
     rfftw_destroy_plan(p);
}

Обратите внимание, что он обрабатывает длины данных, которые не являются четными. Обратите особое внимание, если указана длина данных, FFTW выдаст вам «корзину», соответствующую частоте Найквиста (частота дискретизации, деленная на 2). В противном случае вы не получите его (т. Е. Последний бин находится чуть ниже Найквиста).

A Пример MATLAB похож, но они выбирают длину 1000 (четное число) для примера: 

N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;

В общем, это может зависеть от реализации (ДПФ). Вам следует создать тестовую чистую синусоидальную волну с известной частотой, а затем убедиться в том, что расчет дает то же число.

0 голосов
/ 19 ноября 2010

Частота = скорость / длина волны.

Длина волны - это расстояние между двумя пиками.

...