fft — переполнение стека дискретного преобразования Фурье

Question

fft — переполнение стека дискретного преобразования Фурье

Я пытаюсь написать простые функции 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;
}

1

c++complextype fft

Решение

Другие решения

Есть библиотека 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

-2

Источник

Accepted Answer

Ваше обратное преобразование Фурье явно нарушено: вы игнорируете аргументы комплексных чисел 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;
}

1