Algorithm to zero pad data before FFT

Question

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

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: 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
#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;
}

Answer 1

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.

Answer 2

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
xpad = [x, zeros(1,(K-1)*length(x))];

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

% Plot padded FFT vs unpadded FFT
figure(1);
plotFFT(xpad,length(xpad),Fs,'b'); hold on;
plotFFT(x,length(x),Fs,'-ro'); hold off;
title('Padded vs. Unpadded FFT'); legend('Padded FFT','Unpadded FFT');
xlabel('Frequency (Hz)'); ylabel('abs( X(f) )');