Как рассчитать частоту с быстрым преобразованием Фурье - PullRequest
Новая
       15

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

0 голосов
/ 16 сентября 2018

Я работаю над проектом датчика пульса. Я получаю данные от датчика следующим образом: voltage without artifacts

Теперь я хотел бы использовать этот код: 

  #include <complex>
#include <iostream>
#include <valarray>

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x)
{
    const size_t N = x.size();
    if (N <= 1) return;

    // divide
    CArray even = x[std::slice(0, N/2, 2)];
    CArray  odd = x[std::slice(1, N/2, 2)];

    // conquer
    fft(even);
    fft(odd);

    // combine
    for (size_t k = 0; k < N/2; ++k)
    {
        Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
        x[k    ] = even[k] + t;
        x[k+N/2] = even[k] - t;
    }
}

// Cooley-Tukey FFT (in-place, breadth-first, decimation-in-frequency)
// Better optimized but less intuitive
// !!! Warning : in some cases this code make result different from not optimased version above (need to fix bug)
// The bug is now fixed @2017/05/30 
void fft(CArray &x)
{
    // DFT
    unsigned int N = x.size(), k = N, n;
    double thetaT = 3.14159265358979323846264338328L / N;
    Complex phiT = Complex(cos(thetaT), -sin(thetaT)), T;
    while (k > 1)
    {
        n = k;
        k >>= 1;
        phiT = phiT * phiT;
        T = 1.0L;
        for (unsigned int l = 0; l < k; l++)
        {
            for (unsigned int a = l; a < N; a += n)
            {
                unsigned int b = a + k;
                Complex t = x[a] - x[b];
                x[a] += x[b];
                x[b] = t * T;
            }
            T *= phiT;
        }
    }
    // Decimate
    unsigned int m = (unsigned int)log2(N);
    for (unsigned int a = 0; a < N; a++)
    {
        unsigned int b = a;
        // Reverse bits
        b = (((b & 0xaaaaaaaa) >> 1) | ((b & 0x55555555) << 1));
        b = (((b & 0xcccccccc) >> 2) | ((b & 0x33333333) << 2));
        b = (((b & 0xf0f0f0f0) >> 4) | ((b & 0x0f0f0f0f) << 4));
        b = (((b & 0xff00ff00) >> 8) | ((b & 0x00ff00ff) << 8));
        b = ((b >> 16) | (b << 16)) >> (32 - m);
        if (b > a)
        {
            Complex t = x[a];
            x[a] = x[b];
            x[b] = t;
        }
    }
    //// Normalize (This section make it not working correctly)
    //Complex f = 1.0 / sqrt(N);
    //for (unsigned int i = 0; i < N; i++)
    //  x[i] *= f;
}


int main()
{
    const Complex test[] = { 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0 };
    CArray data(test, 8);

    // forward fft
    fft(data);

    std::cout << "fft" << std::endl;
    for (int i = 0; i < 8; ++i)
    {
        std::cout << data[i] << std::endl;
    }

    }
    return 0;
}

Чтобы вычислить частоту моего сердцебиения. Проблема, которая у меня есть сейчас, заключается в том, чтобы понять результат, который я получаю: быстрое преобразование Фурье (4,0) (1, -2,41421) (0,0) (1, -0,414214) (0,0) (1,0.414214) (0,0) (1,2.41421)

Что это? Я думаю, что это амплитуда и фаза, но я не знаю, как рассчитать из этой частоты.

Частота импульса составляет от 0,33 Гц до 3 Гц. Разрешение этого кода достаточно хорошее? Мне никогда не приходилось работать с преобразованием Фурье

Спасибо за помощь. Я с нетерпением жду ваших ответов.

1 Ответ

0 голосов
/ 17 сентября 2018

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

Частота определяется положением в массиве. Элемент с индексом i предназначен для частоты дискретизации, умноженной на i / N , где N - количество элементов в БПФ .

...