Из ваших вычислений time(u) Я бы сказал, что ваш период выборки Ts равен 2*pi/1000 [s], что приводит к Fs = 1/Ts = 1000/(2*pi) [Hz]. Аналоговая частота f0 вычисляемого вами синуса будет
1*t = 2*pi*f0*t [radians]
f0 = 1/(2*pi) [Hz]
Обратите внимание, что Fs >> f0.
В цифровом домене частота всегда охватывает 2*pi [radians] (это может быть [-pi,pi) или [0,2*pi), НО Eigen вернуть последнее). Таким образом, вам нужно последовательно разделить диапазон [0,2*pi) на N бинов. Например, если индекс равен k, соответствующая нормализованная частота равна f=2*pi*k/N [radians].
Чтобы узнать, какая аналоговая частота f соответствует каждому нормированному интервалу частот, вычислите f = (fs*k/N) [Hz], где fs - частота дискретизации.
О масштабировании и полном спектре из Eigen FFT do c:
1) Масштабирование: другие библиотеки (FFTW, IMKL, KISSFFT) не выполняют масштабирование, поэтому происходит постоянное усиление после прямого и обратного преобразования IFFT (FFT (x)) = Kx; это сделано, чтобы избежать умножения вектора на значение. Недостатком является то, что алгоритмы, которые работали правильно в Matlab / octave, ведут себя не так, как это было реализовано в C ++. Чем отличается Eigen / FFT: выполняется обратимое масштабирование, поэтому IFFT (FFT (x)) = x.
2) Полуспектр реального БПФ: Другие библиотеки используют только половину частотного спектра (плюс один дополнительный образец для Бункера Найквиста) для реального БПФ, другая половина - сопряженная симметрия c первой половины. Это экономит им копию и немного памяти. Недостатком является то, что вызывающий должен иметь специальный лог c для количества бинов в комплексе против реального. Чем отличается Eigen / FFT: полный спектр возвращается из прямого преобразования. Это облегчает общее шаблонное программирование c, устраняя отдельные специализации для реального против сложного. При обратном преобразовании фактически используется только половина спектра, если тип выходного сигнала является действительным.
Итак, вы должны ожидать усиления, просто сделайте тест ifft(fft(x)) == x (проверяется как "мощность ошибки") << «Мощность сигнала»). Вы можете разделить на <code>N, чтобы получить нормализованную версию.
С другой стороны, два всплеска, которые вы видите, являются следствием пункта 2. Графики, которые вы публикуете выше, являются только одной стороной преобразования, другая сторона симметрична c, если сигнал реальный. Вы можете отбросить верхнюю половину вывода.
Этот код:
#include <eigen/Eigen/Dense>
#include <eigen/unsupported/Eigen/FFT>
#include <complex>
#include <cmath>
#include <iostream>
#include <fstream>
unsigned const N = 1000; //
double const Fs = 32; // [Hz]
double const Ts = 1./Fs; // [s]
const double f0 = 5; // [Hz]
std::complex<double> f(std::complex<double> const & t){
return std::sin(2*M_PI*f0*t);
}
int main(){
std::ofstream xrec("xrec.txt");
Eigen::VectorXcd time(N);
Eigen::VectorXcd f_values(N);
Eigen::VectorXd freq(N);
for(int u = 0; u < N; ++u){
time(u) = u * Ts;
f_values(u) = f(time(u));
freq(u) = Fs * u / double(N);
}
Eigen::FFT<double> fft;
Eigen::VectorXcd f_freq(N);
fft.fwd(f_freq, f_values);
for(int u = 0; u < N; ++u){
xrec << freq(u) << " " << std::abs(f_freq(u)) << "\n";
}
}
генерирует xrec.txt. Затем вы можете использовать этот скрипт gnuplot для генерации фигуры:
set key off
set grid
set output "figure.png"
set xlabel "Frequency [Hz]"
plot [-1:34] [-10:500] "xrec.txt" with impulses, "xrec.txt" with points pt 4
На рисунке вы можете увидеть два всплеска с частотой 5 и 27 Гц, как и ожидалось из этого кода. Я изменил значения, чтобы лучше видеть, что происходит, просто попробуйте другие.
В стиле графиков, которые вы показываете, x- диапазон осей [0,16) вместо [0,32), как на этом графике, но, поскольку ваш сигнал реален, спектр симметричен c, и вы можете отбросить половину его.