Why to pad zeros at the middle of sequence instead at the end of the sequence?

Question

One of the implementation of Bluestein algorithm as show in the below link

Bluestein's Algorithm

[conj(W), zeros(1, L-2*N+1), conj(W(N:-1:2)) ]

padding the zeros at the middle of the sequence, eventually gave the correct results.

But padding zeros just at the end does not give the correct results.

Question: When to pad the zeros at the middle of the sequence and when at the end?

Answer 1

When working with the DFT and IFFT we can zero pad the signal, which serves to interpolate new samples in the other domain. We will often see this applied with padding at the end of the sequence or alternatively in the middle of the sequence as referenced in the application linked by the OP. Below I answer the question as to why we would consider padding in the center of the sequence and the implications of doing that. Specific to how that is used in an implementation of Bluestein’s algorithm I have detailed in another post here.

In general when working with the Discrete Fourier Transform (and inverse) , if we don't pad in the middle of a sequence, a linear phase will be introduced in the resulting transform. In applications where our only concern is with the interpolated magnitude of the result, this would be of no consequence.

This answer explains the general considerations and motivations for when we would want to insert zeros in the middle of a sequence in either the time or frequency domain rather than zero padding at the end. Padding a sequence with zeros in one domain in either location (middle or end), interpolates samples in the other domain. This interpolation can be accomplished without introducing additional phase distortion in the other domain by padding in the proper "middle" of the sequence rather than at the end. In many cases this phase distortion is inconsequential since it is a linear phase: A linear phase in the frequency domain is a time shift or delay in the time domain. Similarly, a linear phase in the time domain is a frequency shift or translation in the frequency domain. Alternatively we can zero pad at the end of the sequence and then correct for the linear phase error in the result which may be more convenient that the approaches outlined here.

Proper "Padding in the Center"

Proper symmetry must be maintained when padding in the center, such that we maintain the same number of "positive" and "negative" domain samples.

Padding in the "true center" for an odd sequence is done by placing the zeros in between $N/2+1$ samples at the beginning and then the remaining $N/2$ samples at the end, as in:

$$[x_0, x_1, x_2, x_3, x_4]$$

$$[x_0, x_1, x_2, 0, 0, x_3, x_4]$$

As shown in the link jomega shared in the comments, the location of the zero padding is clear when we consider the alternate and equivalent positive and negative indexing as:

Values: $[x_0, x_1, x_2, x_3, x_4]$

Indexes: $[0, 1, 2, 3, 4]$

Is the same as the following given the periodicity property of the DFT:

Values: $[x_0, x_1, x_2, x_3, x_4]$

Indexes: $[0, 1, 2, -2, -1]$

Thus a zero insert after index 2 can increase the time duration in both the positive and negative direction which serves to not introduce any additional delay:

Values: $[x_0, x_1, x_2, 0, 0, x_3, x_4]$

Indexes: $[0, 1, 2, 3, 4, -4, -3, -2, -1]$

For even sequences, the center bin is shared between the "positive" and "negative" domain, and therefore must be split in complex conjugate halves if not zero. To pad in the center for even sequences, we must split the bin located at $n=N/2$ (for $n=0\ldots N-1$) into complex conjugate halves. For example, with $N=5$, the sample $x_3$ is the shared sample that is right on the boundary between what would be considered the positive domain samples and negative domain samples, and the zero padding would be done as follows:

$$[x_0, x_1, x_2, x_3, x_4, x_5]$$

$$[x_0, x_1, x_2, x_3/2, 0, 0, (x_3/2)^*, x_4, x_5]$$

Where $(x_3/2)^*$ represents the complex conjugate of $x_3/2$.

Related questions with additional details related to the proper splitting are here and here.

Intuition for Padding in the Time Domain

Due to the reciprocity in the DFT, similar considerations in the frequency domain would apply in the time domain by swapping maximum time with maximum frequency. The more detailed frequency domain explanation is given to be more intuitive for anyone familiar with sampling. However, in general for either domain, padding in the center of a sequence will increase the "length" in either domain without distorting the original samples. ("Length" implying the sampling frequency in the frequency domain, or the time duration in the time domain). Below is a simple time domain example, followed by a more detailed frequency domain explanation where we see the same property holds and why.

Consider the time domain sequence given by:

$$x[n] = [1, -1, 1, 1, -1]$$

The DFT of this sequence is:

$$X[k] = [ 1, -1.236, 3.236, 3.236, -1.236]$$

The Discrete Time Fourier Transform (which the DFT samples lie on) is plotted below together with the selection given as $X[k]$ above.

Note that the result $X[k]$ for the particular time domain sequence used is real (due to the symmetry I chose in the sequence and that the samples are real). Likewise the DTFT plotted above is completely real.

Only when we pad zeros in the proper center of the sequence, will the result continue to be real and fall exactly on the plot above. Padding zeros to the end or off-center will result in the samples from the DTFT having the same magnitude as above but with introduced phase offsets (so introduces a linear phase distortion; linear as the phase introduced is proportional to the frequency index).

Intuition for Padding in the Frequency Domain

The N-sample DFT typically has the bins extending from bin 0 to bin N-1 corresponding to the frequencies of DC up to nearly the sampling rate (1 bin less than what would correspond to the sampling rate. Due to periodicity in the DFT, these are equivalent to the DFT bins corresponding to the frequencies in the first Nyquist zone ($|f|< f_s/2$ where $f_s$ is the sampling rate) and as mapped with the fftshift function in MATLAB/Octave and Python scipy.signal.

If we wish to increase the sampling rate associated with the DFT samples while maintaining all signals in the first Nyquist zone, we should then pad the frequency domain DFT result in the center.

This may be made clearer by observing the spectrum plots below showing the representation of a continuous-time (CT) sinusoid in the DFT, along with the spectrum of the sampling process and the resulting spectrum of the discrete-time (DT) sinusoid. If these graphics aren't completely clear, I add more detail on their background at the end of this post. We see how the DFT covers the frequency range from DC to the sampling rate, but also the periodicity in the DFT result such that if we moved the shaded region to the left or right, we could still convey all the information contained about the original signal: the signal component at the upper end of the DFT equivalently represents the "negative frequency" components in the signal.

If we only wanted to indirectly create a new DFT that would be equivalent to sampling the original signal at a higher sampling frequency, then we should zero pad in the middle of the DFT as demonstrated in the graphic below.

Side note: in the actual DFT outputs for the DFT of a sinusoidal tone we would also see many additional non-zero samples as "spectral leakage" except for the convenient case that the sampling rate is an integer multiple of the frequency of the sinusoidal tone.

If the above plots are confusing, please see this post for more background information. Basically the sampling of a signal is the process of multiplying a signal with periodic impulses in time. The Fourier Transform of periodic impulses in time is periodic impulses in frequency (which is what we see in the graphics above). Multiplication in the time domain is identical to convolution in frequency domain. The spectrum for the discrete-time (DT) sinusoid is the result of convolving the CT Sinusoid with the spectrum of the sampling process.

For further examples of this specific to the Bluestein algorithm in this question and the DFT used for efficient circular convolution, and how the zero padding in the middle of a sequence can be done, please see this other post.

Answer 2

The phase zero reference of most all FFT implementations is the first element of the input vector.

If you want the data to remain continuous at and around the phase zero reference, and (circularly) symmetric around the phase zero reference (so that the even portion of the input is represented in the real or cosine result component, and the odd portion is represented in the imaginary or sine result component), then the padding needs to be (circularly) opposite the phase zero reference. Thus around the middle of the input vector, which is (circularly) opposite the first element. Then this padding won't distort the phase (or the even/odd ratio of the input). As it would if the padding were not symmetrically opposite, but unsymmetrically to one side or the other of the phase zero reference.

All the circularly stuff is because all the basis vectors of the FFT are circularly symmetric. But starts with a phase of zero at the first element (and wraps both ways, forward and backwards around the vector, circularly), in the canonical FFT and IFFT implementations.

If you are doing a DFT interpolation, you might not want to corrupt the phase. (e.g. if you want a real and even input to produce a strictly real result, etc.)

Answer 3

I think the other answers unnecessarily obfuscate the matter.

There's nothing special about center padding for sake of convolution. We can pad at any other location, and the unpadded result will be exactly the same - all that changes is the unpadding indices.

Dan answers the more general question of the purpose of center-padding the DFT, which is DFT-upsampling. An alternate perspective is that it's the inverse of subsampling. Regardless that's immaterial to convolution, which is the question's context.

The algorithm

Taking another look, here's what it does

1. Generate W

2. A = [conj(W), 0, flip(conj(W))]

3. B = x * W

4. Right-pad both A and B to length L (done by fft(, L)). Note, x is length N, and so is A (or just suppose it is, I can't confirm because the page won't load for some reason)

5. out = ifft(fft(A_rightpad) * fft(B_rightpad))

6. out = out[:N]. Unpad is 0 to N because we convolved over x * W which was from 0 to N. The "center padding" simply centers the convolving kernel at 0, such that

   • out[0] corresponds to sum(B * A_centered_at_0)
   • out[1] corresponds to sum(B * A_centered_at_1)
   • ...

That's all. It's not about phase or interpolating or upsampling, especially if the example is x0 x1 x2 x3 -> x0 x1 0 0 x2 x3, which changes the very input. Suppose the convolving kernel is two samples long, [h0 h1], and we look at direct convolution where no padding is needed:

conv_at_1_version_1 = x0 * h0 + x1 * h1
conv_at_1_version_2 = x0 * h0 + x1 * h1

conv_at_2_version_1 = x1 * h0 + x2 * h1
conv_at_2_version_2 = x1 * h0 + 0  * h1

Padding must be done outside the sequence, not between the samples.

This is not padding

x0 x1 0 0 x2 x3

just to be crystal clear.

This is not padding, animated

Convolve $x$ with $h$, where

x2 x3 0 0 x1 x2

x0 x1 0 0 x2 x3

Minimal step-by-step

N=4, L=8, len(W)=2

# out[0]
x0 x1 x2 x3 x4 0  0  0  0
w2 w3 0  0  0  0  0  w0 w1

# out[1]
x0 x1 x2 x3 x4 0  0  0  0
w1 w2 w3 0  0  0  0  0 w0

...

which is

out[0] = x0*w2 + x1*w3
out[1] = x0*w1 + x1*w2 + x2*w3
...

so what we unpad in the end, 1 to N, will contain the centerings of the convolution kernel (concat conj W etc) over the input.

Code ("pad location doesn't matter")

import numpy as np
from numpy.fft import fft, ifft

for N_seg in (8, 9, 10):
    x0 = np.random.randn(N_seg) + 1j * np.random.randn(N_seg)
    x1 = np.random.randn(N_seg) + 1j * np.random.randn(N_seg) 
    z  = np.zeros(N_seg)
    y  = np.random.randn(N_seg*4) + 1j * np.random.randn(N_seg*4)
    
    xp0 = np.hstack([x1, z, z, x0])  # zeros centered
    xp1 = np.hstack([z, x0, x1, z])  # signal centered
    xp2 = np.hstack([x0, x1, z, z])  # zeros right
    
    out0 = ifft(fft(xp0) * fft(y))
    out1 = ifft(fft(xp1) * fft(y))
    out2 = ifft(fft(xp2) * fft(y))
    
    assert np.allclose(out0, np.roll(out1, 2*N_seg))
    assert np.allclose(out0, np.roll(out2, 3*N_seg))

Answer 4

close all;
clear all;

x = [2 5 7 9 6];
h = [1 4 7 3 7];

y0 = conv(x, h)

N = length(x);
M = length(h);

L = N + M - 1;
L = 2^nextpow2(L);

xzp = [x, zeros(1,L-N)];
hzp = [h, zeros(1,L-M)];

y1  = ifft( fft(xzp).*fft(hzp) )
y2  = conv(xzp, hzp)

xcp = [x(1:ceil(N/2)), zeros(1,L-N), x(ceil(N/2)+1:N)];
hcp = [h(1:ceil(N/2)), zeros(1,L-N), h(ceil(N/2)+1:N)];

y3  = ifft( fft(xcp).*fft(hcp) )
y4  = conv(xcp, hcp)