# Algorithm to zero pad data before FFT

I've seen some other questions on here about zero padding (like this one) but I'm still a little confused about my situation.

I'm attempting to zero pad my input time series data, so I get an interpolated frequency spectrum. I had initially tried this with NumPy's FFT package, and I checked my algorithm on generated data to see if it works. I generated sine waves of known frequencies, and checked to see what the differences were between the actual, unpadded and padded estimates for frequencies. For frequencies that sat exactly in an FFT bin, the unpadded and padded estimates agreed exactly (as they should).

But I need to implement this in C (using the FFTW library), but the same test fails in the algorithm I implemented. What I'm doing is -

• Read in time series (generated data)
• Perform unpadded FFT, obtain frequency estimate (by looking at the bin with the maximum amplitude)
• To the end of the time series, add zeros (I'm using 10 times as many zeros).
• Take an FFT of the padded array, and obtain frequency estimate.
• Compare the two frequency estimates.

Turns out the two frequency estimates are very different. This is because in my C code, there seem to be multiple places where the value of the FFT matches the maximum.

I think this is because of how I'm padding zeros to the array. Question: How does one correctly pad zeros to an array such that the FFT yields similar results to the unpadded array?

EDIT: Relevant code (Python and 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

npt     = xSize * ySize

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

freqArr = np.fft.fftfreq(ySize) # Frequency bins

datmax  = np.amax(np.abs(dFFT[:,:xSize/2])) # Find magnitude and pos of max val
datwhere = np.where(np.abs(dFFT[:,:xSize/2]) == datmax)

freqD = freqArr[datwhere[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)


C: Note - I've used a few user written functions, but I've tested them and I'm fairly certain that they're not the problem. I'm including them in here so you guys have a working piece of code, though. How I compiled it - 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

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;
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              ii, jj, index;

double           dum;
double           xf, yf, ixf, iyf, xpd, ypd;

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;

npt       = xSize * ySize;
data      = calloc(npt, sizeof(*data));
fftdata   = calloc(npt, sizeof(*fftdata));
freqArr   = calloc(npt, sizeof(*freqArr));

t2f     = fftw_plan_dft_r2c_2d(xSize, ySize, data, fftdata, 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;
}
}
for(ii = ySize; ii < pady; ii++) {
for(jj = xSize; jj < padx; jj++) {
}
}

fftw_execute(t2f);

minmaxArrayC(0, npt, fftdata, &dum, &maxfft);

fftFreqD(ySize, 1, freqArr);

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++) {
}
}

printf("FFT Max: % 5.5lf Where: % 5.5lf\n", maxfft, posmax);

return 0;
}

• Is it possible for you to show us some of the code (say the equivalent bits of Python and C)? Also, why are you doing so much padding? Usually the padding is from the end of your data to the nearest power of the radix of your FFT (e.g 2). Aug 19 '13 at 13:34
• I'll edit in my code shortly! I'm padding for interpolation rather than to get a power of two FFT. I (perhaps naively) assumed that a larger amount of padding will give me better interpolation in the fourier domain. Aug 19 '13 at 13:50
• If you want to interpolate, you'll need to (low pass) filter and downsample - that may be the cause of your difficulties. Aug 19 '13 at 14:12
• @Kitchi You are correct in your assumption that increasing the padding will also increase your resolution in the frequency domain. This result is easier to see when you think in terms of the DFT rather than the FFT. Aug 19 '13 at 16:42
• @TomKealy I don't think that low-pass filtering and downsampling will give the desired result in this case. Upsampling and low-pass filtering will also not be helpful because the peaks of the FFT will not move to their correct locations as a result of these operations. I think that Kitchi has the right idea here, there is just something amiss in the implementation, or in the interpretation of the results. Aug 19 '13 at 16:48

Zero-padding data (that has not been previously non-rectangularly windowed) adds more visible rectangular window artifacts to non-periodic-in-aperture signals. You might be looking at all the local maxima of the ripples of the resulting Sinc function, which are alongside a main lobe maxima.

• I'm hunting for the largest value in the spectrum, and looking at the corresponding bin. If it's the sidelobes that I'm looking at, they won't be as large as the main lobe no? Aug 20 '13 at 4:18
• The two negations are somehow confusing. Zero-padding works only if you have applied any window that fades out the signal to zero. Just cutting out the signal (rectangular window) creates an edge on the continuation point, which will distrub the spectrum. Jun 9 '20 at 19:36

There are instances when some care must be taken when zero padding (ie preservation of signal parity). However, it doesn't seem like it should make a difference in your problem.

There are two things I can think of that might affect your result, or at least your interpretation of the result. The first is scaling; zero padding will affect the average power of your signal. The second is that depending on the number of zeros you choose to add to the end, you can alter the locations of the bin centers.

The MATLAB code below shows an example of what one might expect to see when zero-padding to improve FFT granularity. Notice that both the peak location and peak amplitude become more accurate in the zero-padded version.

Fs = 48e3;                          % Sampling rate (Hz)
f0 = 5e3;                           % Center frequency
N = 2^7;                            % Samples in time series
K = 8;                              % Zero-pad factor

% Generate N time indices
t = (0:N-1)/Fs;

% Create a sinusoid signal of length N*Fs at f0 Hz
x = cos(2*pi*f0*t)/N;

% Create a zero-padded version of x

% Helper function for plotting FFT
plotFFT=@(x,N,Fs,spec) plot(Fs*(0:N-1)/N-Fs/2, abs(fftshift(fft(x))),spec);