Article

Finding the Component Frequencies in a Composite Sine Wave

Use 1D fast Fourier transform to compute the frequency components of a signal.

Overview

Accelerate’s vDSP module provides functions to perform 1D fast Fourier transforms (FFTs) on vectors of data, such as audio signals. The following shows an input signal (at left), and its frequency domain equivalent (at right) after transforming the signal with a forward FFT.

Diagram showing, on the left, a plot of the original signal, and, on the right, a plot of the component frequencies as a series of vertical lines.

You can inspect the frequency-domain data of a forward FFT to compute the individual sine wave components of a composite wave. The technique described in this article is applicable in many digital signal processing applications; for example, finding the dominant frequencies in a dual-tone multi-frequency (DTMF) signal or removing noise from a signal.

In this article you'll learn how to create a composite sine wave, perform a forward FFT, inspect the result to find its components, and perfom an inverse FFT to recreate the original signal.

Create the Composite Signal

Start by creating an array to generate a signal. You’ll need to populate an array of single-precision values using the sums of a series of sine waves. You define the frequencies as the number of cycles per n. The highest measureable frequency, known as the Nyquist frequency, is the element with index n/2, which is 1023 in a zero-based array.

The following creates the array, signal, containing 10 component sine waves:

let n = vDSP_Length(2048)

let frequencies: [Float] = [1, 5, 25, 30, 75, 100,
                            300, 500, 512, 1023]

let tau: Float = .pi * 2
let signal: [Float] = (0 ... n).map { index in
    frequencies.reduce(0) { accumulator, frequency in
        let normalizedIndex = Float(index) / Float(n)
        return accumulator + sin(normalizedIndex * frequency * tau)
    }
}

The following image is a visualization of composite sine waves in signal:

Diagram showing a plot of the orignal signal in the time domain.

Convert Real Data to Interleaved Complex Data

To prepare the real values in signal for use with the FFT functions in vDSP, convert them to an array of complex numbers (represented by DSPComplex structures). Define the real components as the even values in signal, and the imaginary components as the odd values in signal:

let observed: [DSPComplex] = stride(from: 0, to: Int(n), by: 2).map {
    return DSPComplex(real: signal[$0],
                      imag: signal[$0.advanced(by: 1)])
}

Convert Interleaved Complex Data to a Split Complex Vector

The FFT function, vDSP_fft_zrop(_:_:_:_:_:_:_:), accepts split complex vectors for input and output. Use vDSP_ctoz(_:_:_:_:_:) to copy the contents of the interleaved representation of signal to a split complex vector.

let halfN = Int(n / 2)

var forwardInputReal = [Float](repeating: 0, count: halfN)
var forwardInputImag = [Float](repeating: 0, count: halfN)

var forwardInput = DSPSplitComplex(realp: &forwardInputReal,
                                   imagp: &forwardInputImag)

vDSP_ctoz(observed, 2,
          &forwardInput, 1,
          vDSP_Length(halfN))

Because each complex value stores two real values, the length of the DSPSplitComplex structure is half that of signal.

Create the FFT Setup

Create a setup object that contains a precalculated weights array of complex exponentials required to perform the FFT operations. The vaues in the weights array simplify the FFT calculation. Creating this setup object can be expensive, so do it only once; for example, when starting your app. After creating the setup object, you can reuse it later, as needed.

The following code creates a setup object suitable for performing forward and inverse 1D FFTs on a signal containing n elements:

let log2n = vDSP_Length(log2(Float(n)))

guard let fftSetUp = vDSP_create_fftsetup(
    log2n,
    FFTRadix(kFFTRadix2)) else {
        fatalError("Can't create FFT setup.")
}

You can use this setup object for similarly sized smaller FFTs. However, using a weights array built for an FFT that processes a large number of elements can degrade performance for an FFT that processes a significantly smaller number of elements.

When you no longer need the setup object, be sure to deallocate its assigned memory. The following code shows how to use the vDSP_destroy_fftsetup(_:) function to destroy the setup object:

defer {
    vDSP_destroy_fftsetup(fftSetUp)
}

Perform the Forward FFT

Use the vDSP_fft_zrop(_:_:_:_:_:_:_:) function to perform an out-of-place forward 1D FFT on the signal data, creating the frequency domain representation of the signal:

var forwardOutputReal = [Float](repeating: 0, count: halfN)
var forwardOutputImag = [Float](repeating: 0, count: halfN)
var forwardOutput = DSPSplitComplex(realp: &forwardOutputReal,
                                    imagp: &forwardOutputImag)

vDSP_fft_zrop(fftSetUp,
              &forwardInput, 1,
              &forwardOutput, 1,
              log2n,
              FFTDirection(kFFTDirection_Forward))

On return, forwardOutputReal contains the DC components of the forward FFT, and forwardOutputImag contains the Nyquist components of the forward FFT.

Compute Component Frequencies in Frequency Domain

The Nyquist components of the forward FFT contains a series of high-magnitude items, rendered as vertical lines in the graph below:

Diagram showing a plot of the orignal signal in the frequency domain.

These high-magnitude items correspond exaclty to the frequencies you specified in the frequencies array:

let componentFrequencies = forwardOutputImag.enumerated().filter {
    $0.element < -1
}.map {
    return $0.offset
}
        
// Prints "[1, 5, 25, 30, 75, 100, 300, 500, 512, 1023]"
print(componentFrequencies)

Recreate the Original Signal

Use an inverse FFT to recreate a signal in the time-domain from the frequency domain-data returned by the forward FFT:

var inverseOutputReal = [Float](repeating: 0, count: halfN)
var inverseOutputImag = [Float](repeating: 0, count: halfN)

var inverseOutput = DSPSplitComplex(realp: &inverseOutputReal,
                                    imagp: &inverseOutputImag)

vDSP_fft_zrop(fftSetUp,
              &forwardOutput, 1,
              &inverseOutput, 1,
              log2n,
              FFTDirection(kFFTDirection_Inverse))

The process for generating an array of real values from the complex result of the inverse FFT is the opposite of the conversion steps described in Convert Interleaved Complex Data to a Split Complex Vector. Use vDSP_ztoc(_:_:_:_:_:) to copy the contents of the split complex FFT result to an interleaved complex vector:

var recreatedObserved = [DSPComplex](repeating: DSPComplex(real: 0, imag: 0),
                                     count: halfN)

vDSP_ztoc(&inverseOutput, 1,
          &recreatedObserved, 2,
          vDSP_Length(halfN))

Use map to populate the final array of single-precision values. Because the forward transform has a scaling factor of 2 and the inverse transform has a scaling factor of the number of items, divide each result by 2 * n:

let scale = 1 / Float(n * 2)
let recreatedSignal = recreatedObserved.map {
    [$0.real * scale, $0.imag * scale]
}.flatMap {
    return $0
}

On return, recreatedSignal is approximately equal to signal.