Take the 2-minute tour ×

I'm not a DSP expert, but I understand that there are two ways that I can apply a discrete time-domain filter to a discrete time-domain waveform. The first is to convolve them in the time domain, and the second is to take the FFT of both, multiply both complex spectrums, and take IFFT of the result. One key difference in these methods is the second approach is subject to circular convolution.

As an example, if the filter and waveforms are both N points long, the first approach (i.e. convolution) produces a result that is N+N-1 points long, where the first half of this response is the filter filling up and the 2nd half is the filter emptying. To get a steady-state response, the filter needs to have fewer points than the waveform to be filtered.

Continuing this example with the second approach, and assuming the discrete time-domain waveform data is all real (not complex), the FFT of the filter and the waveform both produce FFTs of N points long. Multiplying both spectrums IFFT'ing the result produces a time-domain result also N points long. Here the response where the filter fills up and empties overlap each other in the time domain, and there's no steady state response. This is the effect of circular convolution. To avoid this, typically the filter size would be smaller than the waveform size and both would be zero-padded to allow space for the frequency convolution to expand in time after IFFT of the product of the two spectrums.

My question is, I often see work in the literature from well-established experts/companies where they have a discrete (real) time-domain waveform (N points), they FFT it, multiply it by some filter (also N points), and IFFT the result for subsequent processing. My naive thinking is this result should contain no steady-state response and thus should contain artifacts from the filter filling/emptying that would lead to errors in interpreting the resulting data, but I must be missing something. Under what circumstances can this be a valid approach?

Any insight would be greatly appreciated

share|improve this question
Here is one such example: pcisig.com/specifications/pciexpress/technical_library/… See Figure 20 (blue curve is before and red curve is after filtering FFT spectrum), and Figure 21 is IFFT of red curve in Figure 20. –  ggkmath Oct 2 '10 at 20:01
Please consider a more specific title. I found it necessary to read all 4 paragraphs to get a vague idea of what you are asking, and even then it's not completely clear. –  finnw Jul 9 '11 at 10:56

4 Answers 4

up vote 7 down vote accepted

The basic problem is not about zero padding vs the assumed periodicity, but that Fourier analysis decomposes the signal into sine waves which, at the most basic level, are assumed to be infinite in extent. Both approaches are correct in that the IFFT using the full FFT will return the exact input waveform, and both approaches are incorrect in that using less than the full spectrum can lead to effects at the edges (that usually extend a few wavelengths). The only difference is in the details of what you assume fills in the rest of infinity, not in whether you are making an assumption.

Back to your first paragraph: Usually, in DSP, the biggest problem I run into with FFTs is that they are non-causal, and for this reason I often prefer to stay in the time domain, using, for example, FIR and IIR filters.


In the question statement, the OP correctly points out some of the problems that can arise when using FFTs to filter signals, for example, edge effects, that can be particularly problematic when doing a convolution that is comparable in the length (in the time domain) to the sampled waveform. It's important to note though that not all filtering is done using FFTs, and in the paper cited by the OP, they are not using FFT filters, and the problems that would arise with an FFT filter implementation do not arise using their approach.

Consider, for example, a filter that implements a simple average over 128 sample points, using two different implementations.

FFT: In the FFT/convolution approach one would have a sample of, say, 256, points and convolve this with a wfm that is constant for the first half and goes to zero in the second half. The question here is (even after this system has run a few cycles), what determines the value of the first point of the result? The FFT assumes that the wfm is circular (i.e. infinitely periodic) so either: the first point of the result is determined by the last 127 (i.e. future) samples of the wfm (skipping over the middle of the wfm), or by 127 zeros if you zero-pad. Neither is correct.

FIR: Another approach is to implement the average with an FIR filter. For example, here one could use the average of the values in a 128 register FIFO queue. That is, as each sample point comes in, 1) put it in the queue, 2) dequeue the oldest item, 3) average all of the 128 items remaining in the queue; and this is your result for this sample point. This approach runs continuously, handling one point at a time, and returning the filtered result after each sample, and has none of the problems that occur from the FFT as it's applied to finite sample chunks. Each result is just the average of the current sample and the 127 samples that came before it.

The paper that OP cites takes an approach much more similar to the FIR filter than to the FFT filter (note though that the filter in the paper is more complicated, and the whole paper is basically an analysis of this filter.) See, for example, this free book which describes how to analyze and apply different filters, and note also that the Laplace approach to analysis of the FIR and IIR filters is quite similar what what's found in the cited paper.

share|improve this answer
Thanks tom10, "(that usually extend a few wavelengths)." This is what I observe on the few test cases I've investigated. That is, the middle 1/3'rd of the waveform appears perfectly stable after IFFT of the product of the two spectrums (e.g. freq conv approach). Is there any theory/math that can show the dividing line where the steady-state response becomes perfectly valid in this approach? That is, I'm looking for a way to predict based on theory what region of the filtered time-domain waveform (e.g. the result from IFFT of the product of the filter FFT and the waveform FFT) is stable. –  ggkmath Oct 3 '10 at 19:38
sigh... now I have to get out a book... to quote form Numerical Recipes: "To summarize — we need to pad the data with a number of zeros on one end equal to the maximum positive duration or maximum negative duration of the response function, whichever is larger. (For a symmetric response function of duration M , you will need only M/2 zero pads.) Combining this operation with the padding of the response rk described above, we effectively insulate the data from artifacts of undesired periodicity." see hebb.mit.edu/courses/9.29/2002/readings/c13-1.pdf around Fig. 13.1.14 –  tom10 Oct 3 '10 at 20:06
At least this answers the specific question about convolutions, though not the general consequences of assumed periodicity of FFTs in general. As I recall, there is some theory for that too, but I don't recall it off-hand. –  tom10 Oct 3 '10 at 20:09
Thanks tom10. However, I probably didn't communicate the problem correctly above. Note, the FILTER length EQUALS the waveform length, both N points long (no zero padding in either array). The link above has the filter length < the waveform length. In my case, where both filter and waveform are N points long, I don't know how 0-padding can help... one could add N zero points to both arrays, FFT both arrays, and take the IFFT of the result to end up with 2N-1 time-response, but even here there is no steady-state response (the 1st N pts are filter filling, and last N pts are filter emptying). –  ggkmath Oct 4 '10 at 2:09
As I stated above, my preliminary results show if I have a filter N points long (no zero padding) and a waveform N points long (no zero padding), take the FFT of both, multiply the complex spectrums, and take IFFT, then look at this time-domain result, the middle 1/3rd (approximately) of the time-domain waveform appears perfectly stable (matches expected result exactly). I would have expected the whole time-domain waveform to be the filter filling up, and 100% overlapped with the filter emptying. So it's nice to see some of the array being stable, but I don't know any theory why this occurs. –  ggkmath Oct 4 '10 at 2:14

Here's an example of convolution without zero padding for the DFT (circular convolution) vs linear convolution. This is the convolution of a length M=32 sequence with a length L=128 sequence (using Numpy/Matplotlib):

f = rand(32); g = rand(128)
h1 = convolve(f, g)
h2 = real(ifft(fft(f, 128)*fft(g)))
plot(h1); plot(h2,'r')

alt text The first M-1 points are different, and it's short by M-1 points since it wasn't zero padded. These differences are a problem if you're doing block convolution, but techniques such as overlap and save or overlap and add are used to overcome this problem. Otherwise if you're just computing a one-off filtering operation, the valid result will start at index M-1 and end at index L-1, with a length of L-M+1.

As to the paper cited, I looked at their MATLAB code in appendix A. I think they made a mistake in applying the Hfinal transfer function to the negative frequencies without first conjugating it. Otherwise, you can see in their graphs that the clock jitter is a periodic signal, so using circular convolution is fine for a steady-state analysis.

Edit: Regarding conjugating the transfer function, the PLLs have a real-valued impulse response, and every real-valued signal has a conjugate symmetric spectrum. In the code you can see that they're just using Hfinal[N-i] to get the negative frequencies without taking the conjugate. I've plotted their transfer function from -50 MHz to 50 MHz:

N = 150000                    # number of samples. Need >50k to get a good spectrum. 
res = 100e6/N                 # resolution of single freq point  
f = res * arange(-N/2, N/2)   # set the frequency sweep [-50MHz,50MHz), N points
s = 2j*pi*f                   # set the xfer function to complex radians 

f1 = 22e6       # define 3dB corner frequency for H1 
zeta1 = 0.54    # define peaking for H1 
f2 = 7e6        # define 3dB corner frequency for H2 
zeta2 = 0.54    # define peaking for H2    
f3 = 1.0e6      # define 3dB corner frequency for H3 

# w1 = natural frequency   
w1 = 2*pi*f1/((1 + 2*zeta1**2 + ((1 + 2*zeta1**2)**2 + 1)**0.5)**0.5)  
# H1 transfer function 
H1 = ((2*zeta1*w1*s + w1**2)/(s**2 + 2*zeta1*w1*s + w1**2))            

# w2 = natural frequency 
w2 = 2*pi*f2/((1 + 2*zeta2**2 + ((1 + 2*zeta2**2)**2 + 1)**0.5)**0.5)  
# H2 transfer function  
H2 = ((2*zeta2*w2*s + w2**2)/(s**2 + 2*zeta2*w2*s + w2**2))            

w3 = 2*pi*f3        # w3 = 3dB point for a single pole high pass function. 
H3 = s/(s+w3)       # the H3 xfer function is a high pass

Ht = 2*(H1-H2)*H3   # Final transfer based on the difference functions

subplot(311); plot(f, abs(Ht)); ylabel("abs")
subplot(312); plot(f, real(Ht)); ylabel("real")
subplot(313); plot(f, imag(Ht)); ylabel("imag")

alt text

As you can see, the real component has even symmetry and the imaginary component has odd symmetry. In their code they only calculated the positive frequencies for a loglog plot (reasonable enough). However, for calculating the inverse transform they used the values for the positive frequencies for the negative frequencies by indexing Hfinal[N-i] but forgot to conjugate it.

share|improve this answer
Thanks eryksun, for your nice plot and comments. I agree with the first paragraph. Regarding the second paragraph, you may be right (I'm not a DSP expert) -- would you have any online or other references where I could learn more about how to apply a transfer function to an FFT spectrum? –  ggkmath Nov 19 '10 at 14:09
Also, I'm not sure which graphs you refer to above. The graphs before p.14 use a simple sinusoidal model for illustration only. The graph on p.27 is a snapshot in time, where the start and end times are arbitrarily chosen (note the curve starts and ends on different values). This routine is applied to real analog signals, where nothing is truly periodic. Period jitter is the difference in each clock period from the average clock period. Another way to think about it: a real analog signal always has some random noise on its clock edges, preventing it from being periodic. –  ggkmath Nov 19 '10 at 14:18
Thanks eryksun, with your analysis above, I tend to agree with you that they forgot the conjugate. –  ggkmath Nov 19 '10 at 21:02
Regarding the 33kHz component, it's introduced by spread-spectrum technology, and this is optional to the standard. Thus, the code can't depend on a periodic signal. –  ggkmath Nov 19 '10 at 21:03
Regarding the Laplace xfr function, the standard provides the required xfr function written in Laplace form. I think your approach sounds reasonable though. Wouldn't both end up producing the same result (or, where do errors get introduced by sampling the Laplace xfr function; why to avoid)? –  ggkmath Nov 19 '10 at 21:17

Although there will be artifacts from assuming that a rectangular window of data is periodic at the FFT aperture width, which is one interpretation of what circular convolution does without sufficient zero padding, the differences may or may not be large enough to swamp the data analysis in question.

share|improve this answer
Thanks hotpaw2, it's certainly possible that some data doesn't get effected much, but such a practice seems so prevalent I have a hard time thinking all data is never effected. Typically the literature will use windowing to reduce leakage (e.g. when the samplings are not perfectly periodic in the window of captured data). However, the issue of filter filling/emptying causing a transient leading to additional artifacts is never discussed, which makes me wonder what I'm missing. –  ggkmath Oct 3 '10 at 17:28
The filter can be assumed to be filled by all the periodic images. It may be a common assumption that the tail of the filter response from the previous images is either small enough, or overlaps exactly as if the data really is periodic with the FFT aperture as an integer multiple the period (sometimes it is). –  hotpaw2 Oct 3 '10 at 18:24

I can shed some light to the reason why "windowing" is applied before FFT is applied.

As already pointed out the FFT assumes that we have a infinite signal. When we take a sample over a finite time T this is mathematically the equivalent of multiplying the signal with a rectangular function.

Multiplying in the time domain becomes convolution in the frequency domain. The frequency response of a rectangle is the sync function i.e. sin(x)/x. The x in the numerator is the kicker, because it dies down O(1/N).

If you have frequency components which are exactly multiples of 1/T this does not matter as the sync function is zero in all points except that frequency where it is 1.

However if you have a sine which fall between 2 points you will see the sync function sampled on the frequency point. It lloks like a magnified version of the sync function and the 'ghost' signals caused by the convolution die down with 1/N or 6dB/octave. If you have a signal 60db above the noise floor, you will not see the noise for 1000 frequencies left and right from your main signal, it will be swamped by the "skirts" of the sync function.

If you use a different time window you get a different frequency response, a cosine for example dies down with 1/x^2, there are specialized windows for different measurements. The Hanning window is often used as a general purpose window.

The point is that the rectangular window used when not applying any "windowing function" creates far worse artefacts than a well chosen window. i.e by "distorting" the time samples we get a much better picture in the frequency domain which closer resembles "reality", or rather the "reality" we expect and want to see.

share|improve this answer
Do you mean sinc function and Hamming or Hann window? –  finnw Jul 9 '11 at 11:00

Your Answer


