Many science and engineering applications require an accurate frequency spectrum or Fourier transform of a signal. The Fourier transform of a sequence of samples of a signal is shown in Equation 1.
This equation is written in dimensionless form, where Ù = ùT, T is the sampling interval, and q\[n\] = q(nT), the nth sample of the signal q(t). The equation also assumes that the signal has a finite duration, so that there are only a total of N contiguous nonzero samples. Q(Ù) is a periodic function of the continuous variable Ù, with a period equal to 2ð.
The usual procedure is to evaluate Q(Ù) at a set of N evenly spaced points in the interval Ù = \[0, 2ð\]. This procedure is called a discrete Fourier transform (DFT) and is usually carried out using an algorithm called the fast Fourier transform (FFT). The DFT gives you the Fourier transform at the points Ù = 2ðk/N (k = 0...N −1). In terms of real frequencies, this gives you a resolution of ù = 2ð/(NT) or f = 1/(NT). For many applications, this resolution may be sufficient.
There are some applications, however, where the location of a spectral peak must be determined with great accuracy. In this case, the resolution provided by a DFT may not be adequate. One way to increase the resolution is to simply pad the samples with extra zeros. For example, to double the frequency resolution, you could add N zeros to the end of the q\[n\] sequence, to give a total sequence length of 2N. The DFT will then evaluate Q(Ù) at 2N points, doubling the resolution.
There's a practical limit to how many zeros you can add to increase the resolution. Often, increased resolution is only necessary within a small frequency interval. Here, a simple and straightforward method would be to evaluate Equation 1 directly over the frequency range you're interested in. In what follows, we give an efficient algorithm for doing this.
Using Euler's identity, we can split Equation 1 into an equation for the real part and an equation for the imaginary part:
We efficiently calculate sin(nÙ) and cos(nÙ) recursively, starting only with sinÙ and cosÙ via the Chebyshev polynomials:
Knowing that for Equations 4 and 5, T0 = 1, T1 = cosÙ, U0 = 1, and U1 = 2cosÙ, we can calculate each Tn and Un with the following formulas:
Now we can calculate a subspectrum over a very small frequency range at an arbitrarily high resolution (much higher than the standard FFT), requiring only cosÙ and sinÙ. The program, hrft, is written in C and is executed as follows:
hrft f1 f2 Nf N sps infile outfile
f1 = starting frequency \[Hz\]
f2 = ending frequency \[Hz\]
Nf = number of frequency points to calculate FT
N = number of samples to read from the input file
sps = samples per second
infile = input file
outfile = output file
The figure shows part of the spectrum, produced by an FFT, from 2260 to 2268 Hz, along with the spectrum produced by the high-resolution Fourier transform (HRFT) over that frequency range.
Both methods used the same 2048-point data set. The FFT features a resolution of 4 Hz/bin, while the HRFT, calculated for 41 points, has a resolution of 0.2 Hz/bin. This is evident in the smoothness of the 41-point HRFT plot, as compared to the blocky three-point FFT plot. The data set has a primary frequency component at 2265 Hz, so we can see that the HRFT method correctly identifies this frequency at its peak, while the peak of the FFT is off by 1 Hz. To get as good resolution using the FFT, the data would have to be padded with zeros to a length of 65,536 points.
If you're only interested in calculating the spectrum over a small frequency range, the HRFT may be more efficient than the FFT. The more narrow the frequency range, the more attractive the HRFT becomes over the FFT. This is especially true when the computation must be performed in an embedded system with limited memory resources.