4
$\begingroup$

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?

$\endgroup$
4
  • $\begingroup$ The algo seems to use convolution. If padded conv is the only goal, then padding location only changes the unpadding location and may be done for programmatic convenience -- see "Alternative convolution?" $\endgroup$ Mar 29, 2022 at 19:51
  • $\begingroup$ In the second past it uses fast convolution. $\endgroup$
    – jomegaA
    Mar 30, 2022 at 1:49
  • 1
    $\begingroup$ But I interested in why to zero-pad at the middle of the sequence and how we are sure y(1:N) is the correct sequence instead of any other indices of y. $\endgroup$
    – jomegaA
    Mar 30, 2022 at 1:51
  • $\begingroup$ Added clarification on unpadding $\endgroup$ Mar 31, 2022 at 17:59

4 Answers 4

4
$\begingroup$

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.

DTFT

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.

Spectrum for lower sampling rate case

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.

Spectrum for higher sampling rate case

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.

$\endgroup$
27
  • 1
    $\begingroup$ Padding at the end of a time domain sequence, for example, will introduce a linear phase vs frequency to the prior result in the frequency domain- but padding in the center will leave the original samples unchanged (and interpolate new samples). So this is answering the question stated "Question: When to pad the zeros at the middle of the sequence and when at the end?" in more general terms for all applications not just the one introduced. I'll add a simple example. $\endgroup$ Mar 30, 2022 at 1:07
  • 1
    $\begingroup$ Unsure what you refer to, the entire convolution, not just unpadded portion, is exactly identical, just circularly shifted, unshiftable via np.roll. I'll add an example later $\endgroup$ Mar 30, 2022 at 17:29
  • 1
    $\begingroup$ Why conj(fft(b))? That's not convolution between a and b nor is it what the algorithm does. If you meant fft(conj(b)), then we're just convolving a with c where c=conj(b) and my point stands $\endgroup$ Mar 30, 2022 at 21:06
  • 1
    $\begingroup$ @OverLordGoldDragon oh yes sorry for that distraction- I confused circular correlation w convolution (it is correlation when conjugated)— but regardless my comment holds for that case as well—- for the case of circular correlation using ifft(fft(a) fft(b)) you get very different results if you change where the zero insert is for the reasons I give. $\endgroup$ Mar 30, 2022 at 21:13
  • 1
    $\begingroup$ The example you've shown isn't "padding", it's changing the signal entirely by permuting its parts $\endgroup$ Mar 30, 2022 at 21:29
2
$\begingroup$

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.)

$\endgroup$
4
  • $\begingroup$ Another trivial question, If I want to convolve (fast convolution) the input with a finite impulse response of a system and also keep the phase information in-fact, should I pad zeros at the center for both the input sequence and as well as impulse reponse? $\endgroup$
    – jomegaA
    Mar 30, 2022 at 7:29
  • $\begingroup$ -1 for "if you want convolution results to be correct" - padding location does not affect the unpadded result of convolution. $\endgroup$ Mar 30, 2022 at 16:06
  • $\begingroup$ @OverLordGoldDragon Edited and I believe a good answer now as written no longer deserving of a downvote? $\endgroup$ Mar 31, 2022 at 10:51
  • $\begingroup$ @DanBoschen The answer explains too little for how many upvotes it has, so I'm keeping it $\endgroup$ Mar 31, 2022 at 17:15
2
$\begingroup$

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

enter image description here

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))
$\endgroup$
23
  • $\begingroup$ Attached a small matlab code snippet as my answer and now I am totally confused. $\endgroup$
    – jomegaA
    Mar 30, 2022 at 18:44
  • $\begingroup$ I added a simple example to the beginning of mine- what am I missing? You only padded one of the two sequences, so makes sense it would only be a rotation in this sense (which is the linear phase in the other domain). I think that rotation is the point of the question/answer-- where we pad rotates the result (delay in time of linear phase in freq). Not sure there is any disagreement here. I think that was hotpaws point as well.. $\endgroup$ Mar 30, 2022 at 21:38
  • $\begingroup$ @DanBoschen (cont'd) Let's make it fully explicit: are you saying $[x[0], x[1], 0, 0, x[2], x[3]]$ is padding? where $x[1]$ is the observation/measurement that follows $x[0]$. $\endgroup$ Apr 1, 2022 at 16:23
  • 1
    $\begingroup$ Ok Understood - I am on my phone but will look closely at what you are referring to in my answer. Thanks $\endgroup$ Apr 1, 2022 at 19:10
  • 1
    $\begingroup$ @DanBoschen Apparently I've misread your negative indexing... yes, we seem to agree. But as written your answer reads as "truth in details" as in "devil in details", there's many statements and notation that suggest the fallacy of center padding a sequence like $[x[0], x[1], 0, 0, x[2], x[3]]$, rather than $x[-1], x[-2]$ for the last two. I guess you didn't object to this example earlier because you thought I referred to the array x rather than a discrete-time formula - but that's part of my problem with your notation, I don't know why not just write $[x_0, x_1, 0, 0, x_{-1}, x_{-2}]$ $\endgroup$ Apr 2, 2022 at 17:44
0
$\begingroup$
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)
$\endgroup$
1
  • $\begingroup$ Perhaps xcp and hcp are totally a different sequence? $\endgroup$
    – jomegaA
    Mar 30, 2022 at 19:27

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.