In Bluestein's algorithm, a DFT is alternatively determined by multiplying the time domain sequence with a "chirp", filtering (convolving) that result with the inverse of the chirp, and then multiplying that result with the same chirp. Doing this will provide the same result as a DFT on the original time sequence.
The motivation for doing this is that the radix-2 FFT is very efficient when the sequence is an integer power of 2, but slower otherwise. The approach considered by the OP is to complete the filtering operation (linear convolution) using "fast convolution" techniques that can be done with longer FFTs that are themselves an integer power of 2.
This uses the property of the Fourier Transform that a product in the frequency domain is convolution in the time domain. When using DFT's, this becomes a circular convolution, so must zero pad in order to get the linear convolution needed. For an original length of $N$ samples, the zero padding must be as such that the new length $M > 2N-1$. Further due to the motivation for using this algorithm, we desire that $M$ be an integer power of 2.
That said, the resulting equivalence of the Bluestein's algorithm to the DFT is summarized below with more complete derivation detailed by Sidney Burrus here. What we note from this are three main points:
- We need to compute a linear convolution one way or another to use Bluestein's Algorithm.
- We need only $N$ samples of the linear convolution result (that is otherwise $2N-1$ long).
- The samples we need can be computed efficiently with $2^n$ length FFT's following the zero padding approach given.
To understand past that, we really need to understand the algorithm's operations. I go through the algorithm step by step below and then show how the convolution is accomplished (as one way to implement the algorithm, in the end we need the first $N$ samples for the linear convolution).
Given a DFT as:
$$X(k) = \sum_{n=0}^{N-1} x[n] W_N^{nk} \tag{1}\label{1}$$
Where $W_N^{nk}$ are the so-called "Twiddle Factors" with $W_N^m$ as the "roots of unity" and a convenient short hand for the exponential $e^{-j2\pi nk /N}$. $W_N$ is a phasor on the complex plane with magnitude $1$ and angle $-2\pi/N$, if we raise that to integer exponents from $0$ to $N-1$, we will get uniformly spaced phasors around the unit circle in steps of angle $2\pi/N$, and each of these if raised to the Nth power will equal $1$ and hence are the "roots of unity".
The FFT is an algorithm (with several variants) to computing the DFT result efficiently, and exactly. Bluestein's algorithm is another that will also compute the DFT result, exactly. What we will see below is an alternate approach to computing the solution to \ref{1}, and get the identical result.
As the link given explains, we can alternatively rewrite the DFT using the substitution:
$$nk = \frac{n^2-(k-n)^2+k^2}{2}$$
Confirming this:
$$n^2-(k-n)^2+k^2 = n^2- k^2 + 2nk - n^2 +k^2 = 2nk$$
Beautiful.
So if we
given as copied from the link with the derivation from just substituting $nk$ with $(n^2-(k-n)^2+k^2)$ as:
$$X(k) = \sum_{n=0}^{N-1} x[n] W_N^{n^2/2}W_N^{-(k-n)^2/2}W_N^{k^2/2} \tag{2}\label{2}$$
The above is clearer to see the substitution made, but I will rewriting \ref{2} for clarity:
$$X(k) = W_N^{k^2/2} \sum_{n=0}^{N-1} (x[n] W_N^{n^2/2})W_N^{-(k-n)^2/2}\tag{3}\label{3}$$
Since $W_N$ is a phasor with magnitude $1$ and angle $-2\pi/N$, we can see how if we raised it to the power $n$ with $n$ as an integer count starting at 0, the result would be a spinning phasor rotating at a constant rate: a fixed tone (which is exactly how the DFT works: a correlation with fixed tones). Raising it to $n^2/2$ instead causes the rate of rotation to continuously increase as the count increases, which is a "chirp". I show the real part of such a chirp below using $N=50$:
The red dots are the actual samples for integer $n$, I plotted it together with a blue trace for continuous time so that we can better visualize the "chirp" which is a frequency ramp with time. Since instantaneous frequency is the time derivative of phase , we also see that chirp will be a linear frequency ramp vs time.
Note the repeated use of the general form $W_N^{m^2/2}$ which we can substitute with $h[m]$ such that $h[m] = W_N^{m^2/2}$. This can help us make sense of the core algorithm by rewriting \ref{3} and adding parenthesis to emphasize the resulting three steps:
$$X(k) = h[k] \bigg(\sum_{n=0}^{N-1} (x[n] h[n]) h^{-1}[k-n] \bigg) \tag{4}\label{4}$$
So the above reveals the alternate approach as the high level algorithm, which we ultimately do efficiently using FFT's that are an integer power of 2:
Step 1) Multiply the complete time domain sequence $x[n]$ by a chirp sequence $h[n]$, resulting in $N$ samples of $y[n]$ for $n=0 \ldots N-1$.
$$y[n] = x[n]h[n]$$
Step 2) Convolve (linear convolution!) this result with the inverse chirp. Since the chirp is simply a unit length phasor rotating clockwise on the complex plane (at an increasing rate), the inverse chirp simply rotates counter-clockwise with the same increase in rotation rate.
$$w[k] = y[k] * h^{-1}[k]$$
Note the generalized equation for linear convolution is given as:
$$y[k] * h^{-1}[k] = \sum_{n=-\infty}^{\infty} y[n] h^{-1}[k-n]$$
In our case $n=0\ldots N-1$, resulting in the form for the convolution which also matches \ref{4}:
$$w[k] = y[k] * h^{-1}[k] = \sum_{n=0}^{N-1} y[n] h^{-1}[k-n]$$
Notice with the result here is key to the OP's question: $k$ is indexed from $0$ to $N-1$, $n$ is also indexed from $0$ to $N-1$, as we effectively sweep through all possible values of $n$ and $k$, the result $h^{-1}[m]=h^{-1}[k-n]$ will range from $h^{-1}[-N+1]$ for $k=0, n=N-1$ to $h^{-1}[N-1]$ for $k=N-1, n =0$. Thus we require an index for the generalized $h^{-1}[m]$ that ranges from $m=-N+1$ up to $m=N-1$. We note that the negative indexes for $h[m]$ and $h^{-1}[m]$ are equal to the positive indexes according to $h[-m]= h[m]$:
$$h[-m] = W_N^{(-m)^2/2}= W_N^{m^2/2} = h[m]$$
This explains the time reversed forms used in populating the inverse chirp vector for the convolution.
Step 3: Multiply the result from the linear convolution $w[k]$ with the chirp waveform to get the final DFT result. Note here importantly that the "chirp" in this case is in the frequency domain instead of the time domain. This final chirp is a rotating phasor as we sweep frequency not time! The magnitude of the frequency is always one and the phase versus frequency increases non-linearly as the frequency increases from DC to the sampling rate.
$$X(k) = w[k]h[k]$$
This last step can be interpreted as a "phase correction"; the magnitude result after step 2 will match the DFT result and this step corrects for the phase introduced by the convolution with original waveform multiplied by the chirp. So in cases where we are only interested in the magnitude of the DFT result (such as when computing power spectral densities), we can avoid the additional computations in this third step entirely! I'll refer to that as "Bluestein Lite" in the future.
Given the functional operation of all three steps above to create identically the DFT result for the original $N$ sample sequence, we can accomplish step 2 very efficiently using FFT's (or compute directly using convolution if we desire), but we must take caution in implementing the results we would get for linear convolution, not the circular convolution that the fft approach directly provides:
$$x \circledast y = \text{ifft}\bigg\{\text{fft}\{x\}\text{fft} \{y\}\bigg\}$$
For this we can get the linear convolution result with proper zero padding. And for all the reasons I gave in this other post regarding the sensitivity of where we zero-pad, changing the zero-padding location in the sequence will affect the alignment of the linear convolution result in step 2 with the product result in step 3. (We could pad differently but then require additional shifting elsewhere). It can further be shown at this point how the approach used in the links provided by the OP will match these steps exactly (and not otherwise).
