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 here 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.

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 that contains 2048 elements.

The following creates the array, called signal, that contains 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.

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.FFT(log2n: log2n,
                              radix: .radix2,
                              ofType: DSPSplitComplex.self) 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.

Create the Source and Destination Arrays for the Forward FFT

The FFT operates on complex numbers, that is numbers that contain a real part and an imaginary part. Create two arrays—one for the real parts and one for the imaginary parts—for the input and output to the FFT operation:

let halfN = Int(n / 2)
        
var forwardInputReal = [Float](repeating: 0,
                               count: halfN)
var forwardInputImag = [Float](repeating: 0,
                               count: halfN)
var forwardOutputReal = [Float](repeating: 0,
                                count: halfN)
var forwardOutputImag = [Float](repeating: 0,
                                count: halfN)

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

Perform the Forward FFT

You use DSPSplitComplex structures to pass pointers to the real and imaginary parts arrays to the FFT transform function. Use withUnsafeMutableBufferPointer(_:) to pass the output arrays to

To perform the forward FFT:

  1. Create a DSPSplitComplex structure to store signal represented as complex numbers.

  2. Use convert(interleavedComplexVector:toSplitComplexVector:) to convert the real values in signal to complex numbers. The conversion stores the even values in signal as the real components in forwardInput, and the odd values in signal as the imaginary components in forwardInput.

  3. Create a DSPSplitComplex structure with pointers to forwardOutputReal and forwardOutputImag to receive the FFT result.

  4. Perform the forward FFT.

The following code shows how to perform the forward FFT using the steps described above.

forwardInputReal.withUnsafeMutableBufferPointer { forwardInputRealPtr in
    forwardInputImag.withUnsafeMutableBufferPointer { forwardInputImagPtr in
        forwardOutputReal.withUnsafeMutableBufferPointer { forwardOutputRealPtr in
            forwardOutputImag.withUnsafeMutableBufferPointer { forwardOutputImagPtr in
                
                // 1: Create a `DSPSplitComplex` to contain the signal.
                var forwardInput = DSPSplitComplex(realp: forwardInputRealPtr.baseAddress!,
                                                   imagp: forwardInputImagPtr.baseAddress!)
                
                // 2: Convert the real values in `signal` to complex numbers.
                signal.withUnsafeBytes {
                    vDSP.convert(interleavedComplexVector: [DSPComplex]($0.bindMemory(to: DSPComplex.self)),
                                 toSplitComplexVector: &forwardInput)
                }
                
                // 3: Create a `DSPSplitComplex` to receive the FFT result.
                var forwardOutput = DSPSplitComplex(realp: forwardOutputRealPtr.baseAddress!,
                                                    imagp: forwardOutputImagPtr.baseAddress!)
                
                // 4: Perform the forward FFT.
                fftSetUp.forward(input: forwardInput,
                                 output: &forwardOutput)
            }
        }
    }
}

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 contain 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 exactly 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, using the frequency-domain data returned by the forward FFT.

To perform the forward FFT:

  1. Create the source of the inverse FFT, with pointers to forwardOutputReal and forwardOutputImag.

  2. Create a DSPSplitComplex structure to receive the FFT result.

  3. Perform the inverse FFT.

  4. Return an array of real values from the FFT result. 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:

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

let recreatedSignal: [Float] = forwardOutputReal.withUnsafeMutableBufferPointer { forwardOutputRealPtr in
    forwardOutputImag.withUnsafeMutableBufferPointer { forwardOutputImagPtr in
        inverseOutputReal.withUnsafeMutableBufferPointer { inverseOutputRealPtr in
            inverseOutputImag.withUnsafeMutableBufferPointer { inverseOutputImagPtr in
                
                // 1: Create a `DSPSplitComplex` that contains the frequency domain data.
                let forwardOutput = DSPSplitComplex(realp: forwardOutputRealPtr.baseAddress!,
                                                    imagp: forwardOutputImagPtr.baseAddress!)
                
                // 2: Create a `DSPSplitComplex` structure to receive the FFT result.
                var inverseOutput = DSPSplitComplex(realp: inverseOutputRealPtr.baseAddress!,
                                                    imagp: inverseOutputImagPtr.baseAddress!)
                
                // 3: Perform the inverse FFT.
                fftSetUp.inverse(input: forwardOutput,
                                 output: &inverseOutput)
                
                // 4: Return an array of real values from the FFT result.
                let scale = 1 / Float(n * 2)
                return [Float](fromSplitComplex: inverseOutput,
                               scale: scale,
                               count: Int(n))
            }
        }
    }
}

On return, recreatedSignal is approximately equal to signal.

See Also

Fourier and Cosine Transforms

Using Windowing with Discrete Fourier Transforms

Multiply signal data by window sequence values to reduce spectral leakage.

Signal Extraction from Noise

Use Accelerate’s discrete cosine transform to remove noise from a signal.

Halftone Descreening with 2D Fast Fourier Transform

Reduce or remove periodic artifacts from images.