Bluestein’s Algorithm or Fourier Transform as Linear Chirp FIR Filter. C Implementation Using the Octave GNU Tool

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(N-1).

We need to find spectral estimates for M frequencies:

See the Figure 1-1:

Figure 2-1: 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(M-1) 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+M-1 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 3-1

Figure 3-1: 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 (N-1) complex multiplications

• To calculate a convolution with
  
  requires (N-1) M complex multiplications and (N-1) M complex additions
  

• At the end of the algorithm it is required to multiply M – 1 times by
  
  that is M-1 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.

1. 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:
   
   

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

3. Calculate the product:
   
   

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

5. 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 6-1:

Figure 6-1: 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 built-in 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

Note: I recommend that you look still at the monographs [1] and [2] on this topic.

7. Download the Bluestein_Support.m and Bluestein_float.c/h

You can download the files:

Bluestein_Support.m
Bluestein_float.c/h

with the button:

8. Literature / References

[1] A.Oppenheim, R.Schafer “Discrete-Time Signal Processing“, published by Pearson Education, Inc, publishing as Prentice Hall, 1999

[2] L.Rabiner, B.Gold “Theory and application of digital signal processing“, Prentice-Hall, Ing Englewood Cliffs, New Jercy 1975

[3] Richard E. Blahut “Fast algorithms for digital signal processing“, Addison Wesley Longman Publishing Co, 1985