Как вычислить частоту данных с использованием БПФ?

Question

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

Большое спасибо.

Answer 1

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

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

Answer 2

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

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

Тогда Синусоидальная частота is f0 = fs*n0/N Hertz.

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

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

Answer 3

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

Frequency_of_Peak = Data_Sample_Rate * Bin_number_of_Peak/Length_of_FFT ; 

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

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

Answer 4

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

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

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

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

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

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

4. Возьмите Иффт

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

Answer 5

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

y = sin(2 pi f t) 

С:

• сигнала времени у
• ф центральной частоты
• т времени

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

Итак, чтобы перейти на частоту, вы можете отказаться от части отрицательной частоты. Он расположен после положительной частотной части. Кроме того, первым элементом массива является dc-offset, поэтому частота равна 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() 

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

Answer 6

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

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

#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). В противном случае вы не получите его (т. Е. Последний бит чуть ниже Найквиста).

MATLAB example похож, но они выбирают длину 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; 

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