Дискретное преобразование Фурье C++

Я пытаюсь написать простые функции DFT и IDFT, которые станут моим ядром для будущих проектов. Проблема в том, что IDFT возвращает значение, отличное от входного значения, и я не могу понять, где ошибка. Ниже моего исходного кода:

vector<double> input;
vector<double> result;
vector<complex<double>> output;

double IDFT(int n)
{
    double a = 0;
    double b = 0;
    int N = output.size();
    for(int k = 0; k < N; k++)
    {
        double value = abs(output[k]);
        a+= cos((2 * M_PI * k * n) / N) * value;
        b+= sin((2 * M_PI * k * n) / N) * value;
    }
    complex<double> temp(a, b);
    double result = abs(temp);
    result /= N;
    return result;
}
complex<double> DFT(double in, int k)
{
    double a = 0;
    double b = 0;
    int N = input.size();
    for(int n = 0; n < N; n++)
    {
        a+= cos((2 * M_PI * k * n) / N) * input[n];
        b+= -sin((2 * M_PI * k * n) / N) * input[n];
    }
    complex<double> temp(a, b);
    return temp;
}

int main()
{
    input.push_back(55);
    input.push_back(15);
    input.push_back(86);
    input.push_back(24);
    input.push_back(66);
    input.push_back(245);
    input.push_back(76);

    for(int k = 0; k < input.size(); k++)
    {
        output.push_back(DFT(input[k], k));
        cout << "#" << k << ":\t" << input[k] << " \t>> abs: " << abs(output[k]) << " >> phase: " << arg(output[k]) << endl;
    }
    for(int n = 0; n < output.size(); n++)
    {
        result.push_back(IDFT(n));
        cout << result[n] << endl;
    }
    return 0;
}

c++ complextype fft
person Roman Kulaha    schedule 03.08.2018    source источник
comment
Если это не просто образовательное упражнение, зачем изобретать велосипед?   -  person Ben Jones    schedule 03.08.2018
comment
fftw.org   -  person Ben Jones    schedule 03.08.2018
comment
Возможно, если я решу скомпилировать его для какого-нибудь микроконтроллера, у меня не будет проблем с какими-то дополнительными библиотеками, а также я должен понимать, как это работает на низком уровне.   -  person Roman Kulaha    schedule 03.08.2018
comment
Почитайте про арифметику с плавающей запятой.   -  person zdf    schedule 03.08.2018
comment
Настоятельно рекомендуем sourceforge.net/projects/kissfft — небольшой, автономный, быстрый.   -  person mtrw    schedule 03.08.2018
comment
@ZDF Похоже, это не проблема с плавающей запятой - я бы ожидал близких, но не точных результатов, если бы это было так.   -  person Ben Jones    schedule 03.08.2018
comment
Спасибо, но я хочу понять, где моя ошибка, мне не нужны другие библиотеки.   -  person Roman Kulaha    schedule 03.08.2018
comment
@RomanKulaha Вы уже использовали для этого отладчик? печатать операторы внутри функций? Должно быть довольно легко понять, что пошло не так, если вы понимаете алгоритм...   -  person Ben Jones    schedule 03.08.2018
comment
Вход: 55 15 86 24 66 245 76 Выход: 254,121 64,0808 42,3543 50,0043 50,0043 42,3543 64,080   -  person Roman Kulaha    schedule 03.08.2018
comment
Для будущих проектов определенно лучше использовать существующие вещи. В сети полно подобных проектов для встраиваемых систем. Например. этот github.com/kosme/arduinoFFT. Я использовал его, и он отлично работает. Если у вас есть образовательная цель, вы должны опубликовать ответ в математическом сообществе math.stackexchange.com.   -  person 4xy    schedule 03.08.2018
comment
@BenJones Меня нет рядом с компьютером. Теперь, когда я вижу результат, похоже, вы правы.   -  person zdf    schedule 03.08.2018

Ответы (2)


arrow_upward
3
arrow_downward

Ваше обратное преобразование Фурье явно нарушено: вы игнорируете аргументы комплексных чисел output[k].

Это должно выглядеть так:

double IDFT(size_t n)
{
    const auto ci = std::complex<double>(0, 1);
    std::complex<double> result;
    size_t N = output.size();
    for (size_t k = 0; k < N; k++)
        result += std::exp((1. / N) * 2 * M_PI * k * n * ci) * output[k];
    result /= N;
    return std::abs(result);
}

Изменить.

Если вы хотите явно разделить реальную и мнимую части, вы можете использовать:

double IDFT(size_t n)
{
    double a = 0;
    size_t N = output.size();
    for (size_t k = 0; k < N; k++)
    {
        auto phase = (2 * M_PI * k * n) / N;
        a += cos(phase) * output[k].real() - sin(phase) * output[k].imag();
    }
    a /= N;
    return a;
}
person Evg    schedule 03.08.2018
comment
В случае, если exp сбивает с толку, это эквивалентно замене value на output[k].real() и output[k].imag() в IDFT() OP соответственно. - person Ben Jones; 04.08.2018

arrow_upward
-1
arrow_downward

Для компьютеров с процессорами Intel:

Есть библиотека Intel-IPP. Он предлагает множество функций с очень высокой производительностью. Очень сложно написать что-то быстрее, чем их функции, из-за используемых ими векторных операций. Попробуйте: https://software.intel.com/en-us/intel-ipp

https://software.intel.com/en-us/articles/how-to-use-intel-ipp-s-1d-fourier-transform-functions

person CoralK    schedule 03.08.2018