我在这里看到了一些关于零填充的其他问题（比如这个），但我对我的情况仍然有些困惑。
我正在尝试对我的输入时间序列数据进行零填充，所以我得到了一个插值频谱。我最初用 NumPy 的 FFT 包尝试过这个，然后我检查了生成数据的算法，看看它是否有效。我生成了已知频率的正弦波，并检查了频率的实际、未填充和填充估计之间的差异。对于恰好位于 FFT bin 中的频率，未填充和填充的估计值完全一致（应该如此）。
但是我需要在 C 中实现这个（使用 FFTW 库），但是在我实现的算法中同样的测试失败了。我正在做的是——
- 读入时间序列（生成的数据）
- 执行未填充的 FFT，获得频率估计（通过查看具有最大幅度的 bin）
- 在时间序列的末尾，添加零（我使用了 10 倍的零）。
- 对填充阵列进行 FFT，并获得频率估计。
- 比较两个频率估计。
事实证明这两个频率估计是非常不同的。这是因为在我的 C 代码中，似乎有多个地方 FFT 的值与最大值匹配。
我认为这是因为我如何将零填充到数组中。 问题：如何正确地将零填充到数组中，以使 FFT 产生与未填充数组相似的结果？
编辑：相关代码（Python 和 C）
Python：
#! /usr/bin/env python
import numpy as np
PI = np.pi
ixFreq = 10.0 # This is the bin that will have the max
iyFreq = 10.0 # amplitude in the fourier plane
xSize = 24 # Define the size of the array
ySize = 60
xPeriod = xSize/ixFreq # Calculate what the period needs to be
yPeriod = ySize/iyFreq
xxFreq = 1./xPeriod # This is the *actual* frequency.
yyFreq = 1./yPeriod
padFac = 10 # Padding factor
xPad = padFac * xSize
yPad = padFac * ySize
npt = xSize * ySize
padnpt = xPad * yPad
xy = np.mgrid[0:ySize,0:xSize] # 2D grid to generate the sine wave
data = np.sin(2.0*PI * (xy[0]/yPeriod + xy[1]/xPeriod)) # Sine wave
dFFT = np.fft.fft2(data) # Regular (unpadded) FFT
padFFT = np.fft.fft2(data, s=[yPad,xPad]) # Padded FFT
freqArr = np.fft.fftfreq(ySize) # Frequency bins
freqPad = np.fft.fftfreq(yPad)
datmax = np.amax(np.abs(dFFT[:,:xSize/2])) # Find magnitude and pos of max val
datwhere = np.where(np.abs(dFFT[:,:xSize/2]) == datmax)
padmax = np.amax(np.abs(padFFT[:,:xPad/2]))
padwhere = np.where(np.abs(padFFT[:,:xPad/2]) == padmax)
freqD = freqArr[datwhere[0][0]]
freqP = freqPad[padwhere[0][0]]
print "Actual Freq: % 5.5lf" % yyFreq # Print everything out
print "Data Max: % 5.5lf" % (np.amax(data))
print "FFT Max: % 5.5lf Where: % 5.5lf\n" % (datmax/np.sqrt(npt), freqD)
print "Pad Max: % 5.5lf Where: % 5.5lf\n" % (padmax/np.sqrt(npt), freqP)
C：注意——我使用了一些用户编写的函数，但我已经对它们进行了测试，我相当肯定它们不是问题所在。不过，我将它们包括在这里，所以你们有一段工作代码。我是如何编译它的-gcc -Wall -o "fftwTest" fftwTest.c -lfftw3 -lm
# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# include <complex.h>
# include <fftw3.h>
#define LARGENUMBER 9999999
#define BADDATA -999999
void fftFreqD(int len, double d, double *freqArr)
{
// Calculate given the bin, what the frequency is. Equivalent of
// numpy's fftfreq function
int ii;
freqArr[0] = 0.; //zero bin is always DC component
if (len % 2 == 0){ // If length of even size
for(ii = 1; ii < len/2; ii++)
freqArr[ii] = ii/(d*len);
for(ii = len/2; ii < len; ii++)
freqArr[ii] = -(len-ii)/(d*len);
}
else{ // Else if odd size
for(ii = 1; ii < (floor(len/2)+1); ii++)
freqArr[ii] = ii/(d*len);
for(ii = (floor(len/2)+1); ii < len; ii++)
freqArr[ii] = -(len-ii)/(d*len);
}
return;
}
void minmaxArrayC(int indx1, int indx2, complex *data, double *mindata,
double *maxdata)
{
// Find the minimum and maximum value of a complex array, by taking
// the absolute value
int ii;
*mindata = LARGENUMBER;
*maxdata = BADDATA;
for (ii = indx1; ii < indx2; ii++) {
if (*mindata > cabs(data[ii]))
*mindata = cabs(data[ii]);
if (*maxdata < cabs(data[ii]))
*maxdata = cabs(data[ii]);
}
}
void generate2DSin(double *outarr,double amp,double freqx, double freqy,
double phase,int lenx, int leny)
{
int ii, jj, index=0;
for(ii = 0; ii < lenx; ii++){
for(jj = 0; jj < leny; jj++)
index = ii*leny + jj;
outarr[index] = amp * sin(2.0*M_PI*(freqx*ii+freqy*jj) + phase);
}
return;
}
int main(int argc, char *argv[])
{
int xSize = 24, ySize = 60;
int padx, pady, padfac;
int npt, padnpt;
int ii, jj, index;
double *data, maxfft, maxpad, posmax=0, posmaxpad=0;
double *paddata;
double dum;
double xf, yf, ixf, iyf, xpd, ypd;
double *freqArr, *freqArrPad;
fftw_complex *fftdata, *padfft;
fftw_plan t2f, t2pad;
ixf = 10.0; // The bin that will have max amplitude in the fourier
iyf = 10.0; // plane.
xpd = xSize/ixf; // The period of the wave, in the x and y directions
ypd = ySize/iyf;
xf = 1./xpd; // The actual frequency of the wave
yf = 1./ypd;
padfac = 32;
pady = padfac * ySize;
padx = padfac * xSize;
padnpt = pady * padx;
npt = xSize * ySize;
data = calloc(npt, sizeof(*data));
fftdata = calloc(npt, sizeof(*fftdata));
freqArr = calloc(npt, sizeof(*freqArr));
padfft = calloc(padnpt, sizeof(*padfft));
paddata = calloc(padnpt, sizeof(*paddata));
freqArrPad = calloc(padnpt, sizeof(*freqArrPad));
t2f = fftw_plan_dft_r2c_2d(xSize, ySize, data, fftdata, FFTW_ESTIMATE);
t2pad = fftw_plan_dft_r2c_2d(padx,pady,paddata,padfft,FFTW_ESTIMATE);
generate2DSin(data, 1., xf, yf, 0., xSize, ySize);
for(ii = 0; ii < ySize; ii++) {
for(jj = 0; jj < xSize; jj++) {
index = ii * xSize + jj;
paddata[index] = data[index];
}
}
for(ii = ySize; ii < pady; ii++) {
for(jj = xSize; jj < padx; jj++) {
index = ii*padx + jj;
paddata[index] = 0.;
}
}
fftw_execute(t2f);
fftw_execute(t2pad);
minmaxArrayC(0, npt, fftdata, &dum, &maxfft);
minmaxArrayC(0, padnpt, padfft, &dum, &maxpad);
fftFreqD(ySize, 1, freqArr);
fftFreqD(pady, 1, freqArrPad);
for(ii = 0; ii < ySize; ii++){
for(jj = 0; jj < xSize; jj++) {
index = ii*xSize + jj;
if(cabs(fftdata[index]) == maxfft){
posmax = freqArr[ii];
}
}
}
for(ii = 0; ii < pady; ii++){
for(jj = 0; jj < padx; jj++) {
index = ii*padx + jj;
if(cabs(padfft[index]) == maxpad)
posmaxpad = freqArrPad[ii];
}
}
printf("FFT Max: % 5.5lf Where: % 5.5lf\n", maxfft, posmax);
printf("pad Max: % 5.5lf Where: % 5.5lf\n", maxpad, posmaxpad);
return 0;
}
