Netlib FFTpack实FFT结果与Scipy Rfft差异及C++读取方法咨询

Let's break down the differences between the two outputs and how to correctly interpret the FFT results from Netlib FFTpack.

1. Core Reason: Conjugate Symmetry of Real FFTs #

When you perform an FFT on a real-valued input, the resulting frequency domain data has a key property: conjugate symmetry. This means the k-th frequency component is the complex conjugate of the (n-k)-th component (where n is the FFT size).

Scipy's rfft leverages this symmetry to save memory: it only returns the first n/2 + 1 complex values (from DC to the Nyquist frequency), since the rest can be inferred from symmetry.

Netlib FFTpack's rfftf on the other hand stores the full set of results in the original input array, using a packed format that includes both real and imaginary parts explicitly.

2. Matching Your Results #

Let's map your C++ output to Scipy's correct (un-truncated) rfft results:

Your Netlib FFTpack output array arr (size 8) follows this format for even n:

• arr[0]: DC component (real part, imaginary part = 0) → matches Scipy's first value 55.394
• arr[1]: Real part of frequency bin 1; arr[2]: Imaginary part of frequency bin 1 → corresponds to Scipy's complex value -5.58618 + 1.66889j
• arr[3]: Real part of frequency bin 2; arr[4]: Imaginary part of frequency bin 2 → corresponds to Scipy's complex value 12.316 + 3.11j
• arr[5]: Real part of frequency bin 3; arr[6]: Imaginary part of frequency bin 3 → corresponds to Scipy's complex value 6.98618 - 0.0991108j
• arr[7]: Nyquist frequency component (real part, imaginary part = 0) → matches Scipy's last value 5.174

Your Scipy code used np.real(rfft(arr, n=8)), which discards the imaginary parts of the complex results—this is why you only saw the real components and fewer values. If you run rfft(arr, n=8) without np.real, you'll get the full complex array that aligns with the Netlib data.

3. How to Correctly Read FFT Values in C++ #

For an even FFT size n:

• DC component: arr[0] (real, imaginary = 0)
• For bins 1 to n/2 - 1:
  • Real part: arr[2*k - 1]
  • Imaginary part: arr[2*k]
• Nyquist component: arr[n-1] (real, imaginary = 0)

For your n=8 case, you can extract the complex components like this:

// DC component
float dc_real = arr[0];
float dc_imag = 0.0f;

// Bin 1
float bin1_real = arr[1];
float bin1_imag = arr[2];

// Bin 2
float bin2_real = arr[3];
float bin2_imag = arr[4];

// Bin 3
float bin3_real = arr[5];
float bin3_imag = arr[6];

// Nyquist component
float nyquist_real = arr[7];
float nyquist_imag = 0.0f;

For odd n, the format is slightly different (there's no Nyquist component, and you'll have (n+1)/2 unique bins), but the core idea of alternating real/imaginary parts for non-DC bins still applies.

内容的提问来源于stack exchange，提问作者dnnagy