1
$\begingroup$

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;
}
$\endgroup$
  • 1
    $\begingroup$ 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). $\endgroup$ – Tom Kealy Aug 19 '13 at 13:34
  • 2
    $\begingroup$ 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. $\endgroup$ – Kitchi Aug 19 '13 at 13:50
  • $\begingroup$ If you want to interpolate, you'll need to (low pass) filter and downsample - that may be the cause of your difficulties. $\endgroup$ – Tom Kealy Aug 19 '13 at 14:12
  • $\begingroup$ @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. $\endgroup$ – nispio Aug 19 '13 at 16:42
  • $\begingroup$ @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. $\endgroup$ – nispio Aug 19 '13 at 16:48
1
$\begingroup$

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) )');

Padded vs. Unpadded FFT

$\endgroup$
  • 4
    $\begingroup$ Good overall answer, but "when zero-padding to improve FFT resolution" is not true. To increase resolution you need to take a larger number of samples into the FFT. A better term would be to say it increases frequency granularity. $\endgroup$ – Tarin Ziyaee Aug 19 '13 at 17:22
  • 1
    $\begingroup$ I think what @user4619 is trying to say is that zero padding adds no new information, it just increases the density of bins across the peak. But two close peaks that weren't resolved in the original FFT won't be resolved in the zero-padded FFT. $\endgroup$ – Kitchi Aug 19 '13 at 17:45
  • 1
    $\begingroup$ @user4619 and Kitchi Point taken. However, in the example above, if I add a second tone to my signal at 5.3 kHz I only see one peak in the N-point FFT. However, zero padding will allow me to discern two distinct peaks. What terminology would you use to account for this? Because I am now able to distinguish two distinct signals that I could not previously distinguish, my inclination would be to still say that I experienced an increase in my ability to resolve frequencies, or "resolution." $\endgroup$ – nispio Aug 19 '13 at 18:06
  • 1
    $\begingroup$ I guess I should start a new topic to discuss this, since there is something that is obvious to everyone but me. $\endgroup$ – nispio Aug 19 '13 at 20:14
  • 1
    $\begingroup$ @nispio It would make a good question, feel free to ask it here. About the specific example like I said, you can say that "Two peaks were visible after interpolation". $\endgroup$ – Tarin Ziyaee Aug 19 '13 at 20:20
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Kitchi Aug 20 '13 at 4:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.