C++中用零填充和1/2滤波器仍存在1D FFTW3导数精度问题

Let's break down why your zero-padding is hurting derivative accuracy (even though it’s supposed to help with higher-order derivatives) and walk through actionable fixes you can test:

Possible Causes & Targeted Fixes #

1. Zero-Padding Introduces Unwanted Boundary Discontinuities #

Your original sine wave is perfectly periodic, but tacking zeros on both ends turns it into a non-periodic signal (sine + flat zero regions). FFTW assumes input is periodically extended, so this sharp jump triggers the Gibbs phenomenon—spurious high-frequency components in the spectrum. Since derivative calculations amplify high frequencies (the frequency-domain derivative is X(k) * jω_k, where ω grows with frequency), these artifacts get blown up and wreck your accuracy.

Fixes:

• If you’re simulating an infinite sine wave, use periodic extension instead of zero-padding (repeat the original sine block to reach 3*nX instead of adding zeros). This keeps the signal periodic, eliminating Gibbs artifacts entirely.
• If zero-padding is necessary (e.g., for better frequency resolution), apply a window function to the original sine wave first. A Hanning or Hamming window will smoothly taper the signal’s ends to zero, reducing the sharp discontinuity. Example code snippet:

  // Apply Hanning window to original nX samples
  for (int i = 0; i < nX; i++) {
      double window = 0.5 * (1 - cos(2 * M_PI * i / (nX - 1)));
      padded_signal[i] = original_signal[i] * window;
  }
  // Fill remaining slots with zeros
  memset(padded_signal + nX, 0, sizeof(double) * (2*nX));
  

2. Incorrect Frequency Axis Calculation for Derivatives #

When scaling to 3*nX, you need to ensure your angular frequencies ω_k match FFTW’s output ordering. FFTW returns frequencies from 0 to fs (sampling rate), with the second half representing negative frequencies. If you miscalculate ω_k (e.g., treating all indices as positive frequencies), your derivative calculation will be fundamentally wrong.

Correct ω_k Calculation:
For a signal of length nX3 = 3*nX, compute the angular frequency for index k like this:

double omega;
if (k <= nX3 / 2) {
    omega = 2 * M_PI * k / nX3;  // Positive frequencies
} else {
    omega = 2 * M_PI * (k - nX3) / nX3;  // Negative frequencies
}

Multiply each FFT coefficient X(k) by I * omega (where I is the imaginary unit) to get the derivative’s frequency-domain representation, then inverse FFT back to the time domain.

3. Frequency Leakage from Non-Integer Multiple of Resolution #

If your original sine wave’s frequency f0 is an integer multiple of the original frequency resolution (Δf_old = fs / nX), but not of the new resolution (Δf_new = fs / nX3), FFT will spread the signal’s energy across multiple frequency bins (leakage). This leakage gets amplified in derivative calculations, leading to worse accuracy.

Fix:
Adjust the padded length nX3 so that f0 is an integer multiple of fs / nX3. For example, if f0 = 5 Hz and fs = 100 Hz, nX3 should be a multiple of 20 (since 100/20 = 5). If you can’t adjust the length, use a window function (as mentioned in point 1) to reduce leakage.

4. Suboptimal FFTW Planner Settings #

FFTW’s default FFTW_ESTIMATE planner prioritizes speed over accuracy, especially for non-power-of-two sizes like 3*nX. Switching to a more thorough planner can improve calculation precision.

Fix:
Use FFTW_MEASURE or FFTW_PATIENT when creating your FFT plan. Example:

fftw_plan plan = fftw_plan_dft_1d(nX3, in, out, FFTW_FORWARD, FFTW_MEASURE);

Note that FFTW_PATIENT takes longer to plan but delivers better performance/accuracy for complex sizes.

Debugging Steps to Validate #

• Check the Spectrum: Run FFT on your padded signal and plot the magnitude. Extra high-frequency spikes (beyond the sine’s fundamental) confirm Gibbs artifacts or leakage.
• Verify Derivative Math: Test with a small, known signal (e.g., a sine wave where you can compute the analytical derivative exactly) and compare the FFT-derived result to the analytical one step by step.
• Compare Padding Methods: For the same input, compute the derivative with no padding, zero-padding, and periodic extension—this will isolate whether the padding method is the root issue.

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