Bluestein’s Algorithm or Fourier Transform as Linear Chirp FIR Filter. C Implementation Using the Octave GNU Tool
1. Abbreviation
DSP – Digital Signal Processing
DFT – Discrete Fourier Transform
IDFT – Inverse Discrete Fourier Transform
FFT – Fast Fourier Transform
FIR – Finite Impulse Response
IIR – Infinite Impulse Response
a + jb – complex number where
(a + jb)* = a – jb complex conjugate operation
2. Introduction
In previous articles I have already written about the spectral estimates of the signal using FFT and the Goertzel algorithm. See the articles:

“Fast Fourier Transform (FFT) and C Implementation Using the Octave GNU Tool“

“Goertzel Algorithm and C Implementation Using the Octave GNU Tool“
In this article I want to consider the Bluestein’s algorithm which is inferior to FFT in terms of complexity, but it allows you to calculate the spectrum at arbitrary equidistant points on the unit circle of the z plane and does not require an input array with a size equal to the power of two:
N – can be a prime number. See the example below
Task statement:
Let there be an input sequence of N samples: x(0), x(1),…, x(N1).
We need to find spectral estimates for M frequencies:
See the Figure 11:
Figure 21: The location of M frequencies for calculating the signal spectrum
It is required to calculate the signal spectrum for M frequencies:
Here the description uses the normalized angular frequency ω, which is related to the real angular frequency by the formula:
Bluestein has done the algorithm in the form of a FIR filter with a chirp signal as an impulse response. I described chirp signal already in the article: “Matched Filter Using Octave GNU Tool“:
The frequency of such a signal changes linearly with time. Using this signal you can build an effective matched filter which is used for example in radar systems.
Sampling of the signal (4):
Such a signal has an argument n^2. Moving on to the complex exponent, we get:
3. Bluestein’s Algorithm
Substituting (1) into (2):
Introducing the notation:
Then:
We use the obvious equality:
Substituting (10) into (9):
Let’s introduce a new notation for the input signal:
Then substituting (12) into (11):
It is obvious that (13) is a convolution. See the (14):
Then (13) can be rewritten as:
Thus the Fourier Transform for M frequencies from an input signal with a length of N samples can be calculated using a filter having a pulse response in the form of a chirp signal:
We need only the first M values of X(0),…, X(M1) at the output of the filter (15). Also given that the input signal has a finite number of N input samples, then in this case a finite pulse response (16) of length N+M1 samples is required. That is a FIR filter. Then (16) can be written in a simpler (refined) form:
(17) determines the impulse response of the filter for negative samples (negative time), which complicates the filter implementation. To correct this shortcoming we redefine (17) for positive time:
Then in this case (13) and (15) can be rewritten as:
Note:
Formulas (18) and (19) define the Bluestein’s algorithm as a FIR filter with chirp signal.
See the Figure 31
Figure 31: Bluestein’s Algorithm
4. Estimation of the Bluestein’s Algorithm Complexity
The algorithm is performed on complex numbers so we will make two remarks:

Complex multiplication is 4 real multiplications and 2 real additions:

Complex addition is 2 real additions:
So there are N samples of the input signal x(n) and we need to make M spectral estimates. Then:

To calculate
by the formula (12), it is required to make (N1) complex multiplications 
To calculate a convolution with
requires (N1) M complex multiplications and (N1) M complex additions

At the end of the algorithm it is required to multiply M – 1 times by
that is M1 complex multiplications
Thus Bluestein’s algorithm requires:
On the other hand the calculation of the spectrum for M frequencies by the formula (2) requires less cost:
Obviously the direct use of Bluestein’s algorithm is ineffective. But using the fast convolution, described in the article: “Fast Fourier Transform (FFT) and C Implementation Using the Octave GNU Tool”, allows us to effectively apply Bluestein’s algorithm.
5. Bluestein’s Algorithm using Fast Convolution
Let’s define:
The L choice will allow us to use the FFT algorithm further. Now let’s define the Fast Convolution algorithm for implementing Bluestein’s algorithm.

Add zeros to the impulse response to get an array of size L:
Calculate the L point DFT of the impulse response and store it in the ROM of the microprocessor before the algorithm start:

Add zeros to the input signal to get an array of size L:
Calculate the L point DFT from the input signal using FFT:

Calculate the product:

Next we calculate the inverse Fourier transform using FFT:
Leaving only M samples, and discarding the rest, we get:
y(N – 1), y(N), y(N + 1), …, y(N + M – 2)

Finally calculate the desired spectrum:
Let’s estimate the complexity of the algorithm:

The calculation
requires
Note: the values
are stored in the ROM of the microprocessor before algorithm start

The L point FFT for calculating (28) requires:

Calculation (29):

The L point IFFT for calculating (30) requires:

Calculation (31)
Thus Blustein’s algorithm with fast convolution requires:
Next for simplicity I will compare the algorithms only by the number of multiplications. Now let’s compare the complexity of (23) and (37). Obviously the above algorithm (Bluestein using fast convolution) makes sense to use if the condition (23) > (37):
6. Example of the Bluestein’s Algorithm with the Fast Convolution
Consider the example I used in demo SW: Bluestein_Support.m and Bluestein_float.c/h. Let us have an input signal of length N = 199 and we need to calculate the spectrum at M = 58 points:
See the Figure 61:
Figure 61: Example. The location of M = 58 frequencies for calculating the signal spectrum
Using (24):
Then we can choose: L = 256
Now let’s check the (38):
The Bluestein_Support.m script (Octave Tool) tests the algorithm and generates the help arrays for C Code

To support FFT in C Code the script generates an array of exponents in the file: fft_float_exp_rotation.txt

The Walsh function with a limited length of N=199 samples is taken as the input test signal. The builtin function hadamard () from the Octave Tool was used for the Walsh function calculation. See the array: testSignal. This array is stored in complexWalshTestSignal.txt to support C Code

To calculate (12) the array: inputFactor is used. The values are calculated in this array:
This array is stored in BluesteinInputFactor.txt to support C Code

To calculate (31) the array: outputFactor is used. The values are calculated in this array:
This array is stored in BluesteinOutputFactor.txt to support C Code

To calculate (26) the array: tilde_H is used. This array is stored in BluesteinImpulseResponseSpectrum.txt to support C Code

The result of the calculation of the Bluestein’s algorithm (31) is stored in BluesteinFastConv

Notes

To check the correctness of the obtained result, the spectrum is calculated using the formula (2) too. See the array ClassicalSpectrum

The Bluestein’s algorithm was still executed using the FIR filter function. See the array BluesteinFIR. This variant is not optimal and is given as an additional test

In my article I used the impulse response of Blustein FIR filter (18) from Oppenheim and Schafer [1]. In the monograph Rabiner and Gold [2] a slightly different version of this impulse response is given:
In this case the (31) will be rewritten as:
I have not described this variant here so as not to confuse the reader. However in the Octave script I implemented this option too. See the BluesteinFastConv1

The C implementation is given in Bluestein_float. c/h

flt_complex_nmb_t flt_complex_exp_table[FFT_LOG2_POINTS_NMB] – exponents for the FFT support from the fft_float_exp_rotation.txt

flt_complex_nmb_t BluesteinInputFactors[N_VALUE] – input factors from the BluesteinInputFactor.txt

flt_complex_nmb_t BluesteinOutputFactors[M_VALUE] – output factors from the BluesteinOutputFactor.txt

flt_complex_nmb_t BluesteinImpulseResponseSpectrum[FFT_POINTS_NMB] – Bluestein’s filter impulse response spectrum from the BluesteinImpulseResponseSpectrum.txt

Bluestein_using_fastConvolution( ) – the function is calculated the Bluestein’s algorithm

main( ) – test of the SW:

flt_complex_nmb_t flt_test_signal[FFT_POINTS_NMB] – Walsh test signal from the complexWalshTestSignal.txt

flt_complex_nmb_t flt_output_spectrum[M_VALUE] – the result of the Bluestein’s algorithm will be saved in the array
