C++调用FFTW之注意

原创已于 2022-03-31 09:59:34 修改 · 3.4k 阅读

9 ·

CC 4.0 BY-SA版权

文章标签：

#c++ #python

于 2022-03-31 09:55:20 首次发布

本文介绍了C++调用FFTW库进行傅里叶变换的关键细节。fftw_plan_dft_1d函数中，FFTW_FORWARD表示正变换，FFTW_BACKWARD表示逆变换；FFTW_MEASURE用于精确计算，FFTW_ESTIMATE则追求速度。务必在执行变换前对fftIN赋值，并确保fftIN和fftOut长度相同且初始值为0。

FFTW库函数的fftw_plan_dft_1d中的第三个输入参数FFTW_FORWARD为正傅里叶变换，FTW_BACKWARD为逆傅里叶变换；第四个输入参数FFTW_MEASURE表示准确计算，FFTW_ESTIMATE表示快速估计。

fftw_plan_dft_1d的执行需要在赋值fftIN之前。

fftIn和fftOut的长度需要一致，默认初始化为0。

int main()
{
	int Nfft = 1024, Ns = 600;
	double fc1 = 0e6;
	double fc2 = 1.2e6;
	double Fs = 3.5e6;
	fftw_complex *fftIn = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Nfft);
	fftw_complex *fftOut = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*Nfft);
	fftw_plan plan = fftw_plan_dft_1d(Nfft, fftIn, fftOut, FFTW_FORWARD, FFTW_MEASURE);
	for (int i = 0; i < Ns; i++)
	{
		fftIn[i][0] = cos(2 * PI*fc1*i / Fs)+ cos(2 * PI*fc2*i / Fs); // real part
		fftIn[i][1] = sin(2 * PI*fc1*i / Fs)+ sin(2 * PI*fc2*i / Fs); // imag part
	}
	fftw_execute(plan);
	fftw_destroy_plan(plan);
	double *fftres = new double[Nfft];
	VectorXd freqpoint = VectorXd::LinSpaced(Nfft, -Nfft / 2, Nfft / 2 - 1);
	freqpoint = freqpoint*Fs / Nfft/1e6;
	for (int i = 0; i < Nfft; i++)
	{
		if (i<Nfft/2)
			fftres[i + Nfft / 2] = 10 * log10(sqrt(fftOut[i][0] * fftOut[i][0] + fftOut[i][1] * fftOut[i][1]));
		else
			fftres[i - Nfft / 2] = 10 * log10(sqrt(fftOut[i][0] * fftOut[i][0] + fftOut[i][1] * fftOut[i][1]));
	}
}