How the Cooley-Tukey FFT Algorithm Works | Part 4 - Twiddle Factors
The beauty of the FFT algorithm is that it does the same thing over and over again. It treats every stage of the calculation in exactly the same way. However, this. “one-size-fits-all” approach, although elegant and simple, causes a problem. It misaligns samples and introduces phase distortions during each stage of the algorithm. To overcome this, we need Twiddle Factors, little phase correction factors that push things back into their correct positions before continuing onto the next stage.
How the Cooley-Tukey FFT Algorithm Works | Part 3 - The Inner Butterfly
At the heart of the Cooley-Tukey FFT algorithm lies a butterfly, a simple yet powerful image that captures the recursive nature of how the FFT works. In this article we discover the butterfly’s role in transforming complex signals into their frequency components with efficiency and elegance. Starting with the 2-point DFT, we reveal how the FFT reuses repeated calculations to save time and resources. Using a divide-and-conquer approach, the algorithm breaks signals into smaller groups, processes them through interleaving butterfly diagrams, and reassembles the results step by step.
How the Cooley-Tukey FFT Algorithm Works | Part 2 - Divide & Conquer
The Fast Fourier Transform revolutionized the Discrete Fourier Transform by making it much more efficient. In part 1, we saw that if you run the DFT on a power-of-2 number of samples, the calculations of different groups of samples repeat themselves at different frequencies. By leveraging the repeating patterns of sine and cosine values, the algorithm enables us to calculate the full DFT more efficiently. However, the calculations of certain groups of samples repeat more often than others. In this article, we’re going to explore how the divide-and-conquer method prepares the ground for the next stage of the algorithm by grouping the samples into specially ordered pairs.
How the Cooley-Tukey FFT Algorithm Works | Part 1 - Repeating Calculations
The Fourier Transform is a powerful tool, used in many technologies, from audio processing to wireless communication. However, calculating the FT can be computationally expensive. The Cooley-Tukey Fast Fourier Transform (FFT) algorithm provides a significant speedup. It exploits the repetitive nature of calculations within the Discrete Fourier Transform (DFT), the mathematical foundation of the FT. By recognizing patterns in the DFT calculations and reusing intermediate results, the FFT vastly reduces the number of operations required. In this series of articles, we will look at how the Cooley-Tukey FFT algorithm works.
Learn to Use the Discrete Fourier Transform
Discrete-time sequences arise in many ways: a sequence could be a signal captured by an analog-to-digital converter; a series of measurements; a signal generated by a digital modulator; or simply the coefficients of a digital filter. We may wish to know the frequency spectrum of any of these sequences. The most-used tool to accomplish this is the Discrete Fourier Transform (DFT), which computes the discrete frequency spectrum of a discrete-time sequence. The DFT is easily calculated using software, but applying it successfully can be challenging. This article provides Matlab examples of some techniques you can use to obtain useful DFT’s.
Model a Sigma-Delta DAC Plus RC Filter
Sigma-delta digital-to-analog converters (SD DAC’s) are often used for discrete-time signals with sample rate much higher than their bandwidth. For the simplest case, the DAC output is a single bit, so the only interface hardware required is a standard digital output buffer. Because of the high sample rate relative to signal bandwidth, a very simple DAC reconstruction filter suffices, often just a one-pole RC lowpass. In this article, I present a simple Matlab function that models the combination of a basic SD DAC and one-pole RC filter. This model allows easy evaluation of the overall performance for a given input signal and choice of sample rate, R, and C.
DAC Zero-Order Hold Models
This article provides two simple time-domain models of a DAC’s zero-order hold. These models will allow us to find time and frequency domain approximations of DAC outputs, and simulate analog filtering of those outputs. Developing the models is also a good way to learn about the DAC ZOH function.
Frequency Formula for a Pure Complex Tone in a DTFT
The analytic formula for calculating the frequency of a pure complex tone from the bin values of a rectangularly windowed Discrete Time Fourier Transform (DTFT) is derived. Unlike the corresponding Discrete Fourier Transform (DFT) case, there is no extra degree of freedom and only one solution is possible.
Pentagon Construction Using Complex Numbers
A method for constructing a pentagon using a straight edge and a compass is deduced from the complex values of the Fifth Roots of Unity. Analytic values for the points are also derived.
Interpolator Design: Get the Stopbands Right
In this article, I present a simple approach for designing interpolators that takes the guesswork out of determining the stopbands.
Minimum Shift Keying (MSK) - A Tutorial
Minimum Shift Keying (MSK) is one of the most spectrally efficient modulation schemes available. Due to its constant envelope, it is resilient to non-linear distortion and was therefore chosen as the modulation technique for the GSM cell phone standard.
MSK is a special case of Continuous-Phase Frequency Shift Keying (CPFSK) which is a special case of a general class of modulation schemes known as Continuous-Phase Modulation (CPM). It is worth noting that CPM (and hence CPFSK) is a...
How the Cooley-Tukey FFT Algorithm Works | Part 1 - Repeating Calculations
The Fourier Transform is a powerful tool, used in many technologies, from audio processing to wireless communication. However, calculating the FT can be computationally expensive. The Cooley-Tukey Fast Fourier Transform (FFT) algorithm provides a significant speedup. It exploits the repetitive nature of calculations within the Discrete Fourier Transform (DFT), the mathematical foundation of the FT. By recognizing patterns in the DFT calculations and reusing intermediate results, the FFT vastly reduces the number of operations required. In this series of articles, we will look at how the Cooley-Tukey FFT algorithm works.
Digital Envelope Detection: The Good, the Bad, and the Ugly
Recently I've been thinking about the process of envelope detection. Tutorial information on this topic is readily available but that information is spread out over a number of DSP textbooks and many Internet web sites. The purpose of this blog is to summarize various digital envelope detection methods in one place.
Here I focus on envelope detection as it is applied to an amplitude-fluctuating sinusoidal signal where the positive-amplitude fluctuations (the sinusoid's envelope)...
A Quadrature Signals Tutorial: Complex, But Not Complicated
Introduction Quadrature signals are based on the notion of complex numbers and perhaps no other topic causes more heartache for newcomers to DSP than these numbers and their strange terminology of j operator, complex, imaginary, real, and orthogonal. If you're a little unsure of the physical meaning of complex numbers and the j = √-1 operator, don't feel bad because you're in good company. Why even Karl Gauss, one the world's greatest mathematicians, called the j-operator the "shadow of...
A Fixed-Point Introduction by Example
Introduction
The finite-word representation of fractional numbers is known as fixed-point. Fixed-point is an interpretation of a 2's compliment number usually signed but not limited to sign representation. It extends our finite-word length from a finite set of integers to a finite set of rational real numbers [1]. A fixed-point representation of a number consists of integer and fractional components. The bit length is defined...
Use Matlab Function pwelch to Find Power Spectral Density – or Do It Yourself
In my last post, we saw that finding the spectrum of a signal requires several steps beyond computing the discrete Fourier transform (DFT)[1]. These include windowing the signal, taking the magnitude-squared of the DFT, and computing the vector of frequencies. The Matlab function pwelch [2] performs all these steps, and it also has the option to use DFT averaging to compute the so-called Welch power spectral density estimate [3,4].
In this article, I’ll present some...
Second Order Discrete-Time System Demonstration
Discrete-time systems are remarkable: the time response can be computed from mere difference equations, and the coefficients ai, bi of these equations are also the coefficients of H(z). Here, I try to illustrate this remarkableness by converting a continuous-time second-order system to an approximately equivalent discrete-time system. With a discrete-time model, we can then easily compute the time response to any input. But note that the goal here is as much to...
Interpolation Basics
This article covers interpolation basics, and provides a numerical example of interpolation of a time signal. Figure 1 illustrates what we mean by interpolation. The top plot shows a continuous time signal, and the middle plot shows a sampled version with sample time Ts. The goal of interpolation is to increase the sample rate such that the new (interpolated) sample values are close to the values of the continuous signal at the sample times [1]. For example, if...
Design IIR Filters Using Cascaded Biquads
This article shows how to implement a Butterworth IIR lowpass filter as a cascade of second-order IIR filters, or biquads. We’ll derive how to calculate the coefficients of the biquads and do some examples using a Matlab function biquad_synth provided in the Appendix. Although we’ll be designing Butterworth filters, the approach applies to any all-pole lowpass filter (Chebyshev, Bessel, etc). As we’ll see, the cascaded-biquad design is less sensitive to coefficient...
Digital PLL's -- Part 1
1. IntroductionFigure 1.1 is a block diagram of a digital PLL (DPLL). The purpose of the DPLL is to lock the phase of a numerically controlled oscillator (NCO) to a reference signal. The loop includes a phase detector to compute phase error and a loop filter to set loop dynamic performance. The output of the loop filter controls the frequency and phase of the NCO, driving the phase error to zero.
One application of the DPLL is to recover the timing in a digital...
A Quadrature Signals Tutorial: Complex, But Not Complicated
Introduction Quadrature signals are based on the notion of complex numbers and perhaps no other topic causes more heartache for newcomers to DSP than these numbers and their strange terminology of j operator, complex, imaginary, real, and orthogonal. If you're a little unsure of the physical meaning of complex numbers and the j = √-1 operator, don't feel bad because you're in good company. Why even Karl Gauss, one the world's greatest mathematicians, called the j-operator the "shadow of...
A Fixed-Point Introduction by Example
IntroductionThe finite-word representation of fractional numbers is known as fixed-point. Fixed-point is an interpretation of a 2's compliment number usually signed but not limited to sign representation. It extends our finite-word length from a finite set of integers to a finite set of rational real numbers [1]. A fixed-point representation of a number consists of integer and fractional components. The bit length is defined...
Minimum Shift Keying (MSK) - A Tutorial
Minimum Shift Keying (MSK) is one of the most spectrally efficient modulation schemes available. Due to its constant envelope, it is resilient to non-linear distortion and was therefore chosen as the modulation technique for the GSM cell phone standard.
MSK is a special case of Continuous-Phase Frequency Shift Keying (CPFSK) which is a special case of a general class of modulation schemes known as Continuous-Phase Modulation (CPM). It is worth noting that CPM (and hence CPFSK) is a...
Digital Envelope Detection: The Good, the Bad, and the Ugly
Recently I've been thinking about the process of envelope detection. Tutorial information on this topic is readily available but that information is spread out over a number of DSP textbooks and many Internet web sites. The purpose of this blog is to summarize various digital envelope detection methods in one place.
Here I focus on envelope detection as it is applied to an amplitude-fluctuating sinusoidal signal where the positive-amplitude fluctuations (the sinusoid's envelope)...
Design IIR Butterworth Filters Using 12 Lines of Code
While there are plenty of canned functions to design Butterworth IIR filters [1], it’s instructive and not that complicated to design them from scratch. You can do it in 12 lines of Matlab code. In this article, we’ll create a Matlab function butter_synth.m to design lowpass Butterworth filters of any order. Here is an example function call for a 5th order filter:
N= 5 % Filter order fc= 10; % Hz cutoff freq fs= 100; % Hz sample freq [b,a]=...Use Matlab Function pwelch to Find Power Spectral Density – or Do It Yourself
In my last post, we saw that finding the spectrum of a signal requires several steps beyond computing the discrete Fourier transform (DFT)[1]. These include windowing the signal, taking the magnitude-squared of the DFT, and computing the vector of frequencies. The Matlab function pwelch [2] performs all these steps, and it also has the option to use DFT averaging to compute the so-called Welch power spectral density estimate [3,4].
In this article, I’ll present some...
The Exponential Nature of the Complex Unit Circle
IntroductionThis is an article to hopefully give an understanding to Euler's magnificent equation:
$$ e^{i\theta} = cos( \theta ) + i \cdot sin( \theta ) $$
This equation is usually proved using the Taylor series expansion for the given functions, but this approach fails to give an understanding to the equation and the ramification for the behavior of complex numbers. Instead an intuitive approach is taken that culminates in a graphical understanding of the equation.
Complex...Design IIR Filters Using Cascaded Biquads
This article shows how to implement a Butterworth IIR lowpass filter as a cascade of second-order IIR filters, or biquads. We’ll derive how to calculate the coefficients of the biquads and do some examples using a Matlab function biquad_synth provided in the Appendix. Although we’ll be designing Butterworth filters, the approach applies to any all-pole lowpass filter (Chebyshev, Bessel, etc). As we’ll see, the cascaded-biquad design is less sensitive to coefficient...
Design IIR Bandpass Filters
In this post, I present a method to design Butterworth IIR bandpass filters. My previous post [1] covered lowpass IIR filter design, and provided a Matlab function to design them. Here, we’ll do the same thing for IIR bandpass filters, with a Matlab function bp_synth.m. Here is an example function call for a bandpass filter based on a 3rd order lowpass prototype:
N= 3; % order of prototype LPF fcenter= 22.5; % Hz center frequency, Hz bw= 5; ...How to Find a Fast Floating-Point atan2 Approximation
Context Over a short period of time, I came across nearly identical approximations of the two parameter arctangent function, atan2, developed by different companies, in different countries, and even in different decades. Fascinated with how the coefficients used in these approximations were derived, I set out to find them. This atan2 implementation is based around a rational approximation of arctangent on the domain -1 to 1:
$$ atan(z) \approx \dfrac{z}{1.0 +...
