Fourier Transform and frequency domain
Digital signals like sound are usually recorded and processed as a time domain stream – sequence of amplitude values in regular time intervals. This representation is easy to record, convert and generally deal with. In case of sound however, the “waveform” gives us almost no information on how it actually sounds. Human ears analyze sound by separating it into individual sine waves and measuring amplitudes of the sine waves (and how they change over time). Can similar process be applied on digital time-domain signal? Yes, it can…
The aim of this article is to show you what the frequency domain is, and how to make sense of frequency-domain-represented data.
Let’s start with the most basic possible situation. Let’s say you have a waveform oscillator with constant frequency. Fourier theory says, that Every harmonic oscillator may be divided into several (sometimes infinite number of) sine wave oscillators, each with distinctive phase, magnitude and frequency which is integer multiple of the base frequency. In music theory these frequencies are called harmonics (first harmonic being the fundamental frequency = pitch) or overtones (+ fundamental frequency… in fact first overtone = second harmonic).
So how do we measure the phase and magnitude of these harmonics? Let’s that with the magnitude first… First thing to do is to cut exactly one period of the signal – this will give as a finite array of values to work with. You may remember from math in high school that when you multiply two sinuses with same frequency you get a DC offset and sinus with double the frequency:
One of the sine will be the “testing” sine (which is simple sin(x) function that peaks at +1) and one is the period of signal. When you average the multiplication product of the two, obviously the sin(2x) goes away (because it has same amount above and below zero) and only the DC component remains ( the 0.5) which can simply by multiplied by 2 to get magnitude of the original signal. (Click the thumbnail to enlarge the image)
This works only if the sine waves are the same frequency. If the sine in the signal is of different frequency, the result will always be zero:
And the coolest part – This works even if the signal consists of multiple sines stacked on top of each other. You can still measure only the amplitude of desired harmonic:
There is a twist though… The signal and testing sine wave must be phase aligned. If the signal is shifted to the left or right, you will not get correct magnitude values – they will range anywhere between positive and negative true value of the magnitude (in fact, it may even be zero if the sines are +- pi/2 phase shifted). That’s quite a bugger… One stupid way to solve this would be to do this many times with different phase shifted test sines and pick the one with highest value (you will also then know the phase of the harmonic). However, there is easier way to do this – do two phase shifted sines and “extrapolate” the phase and magnitude from the two measured values – in fact most reasonable is to pick sine and cosine as test waves. Also notice, that when measured value is max or min for sine it is zero for cosine and vice versa.
You can plot these two values into Cartesian graph as a vector. The length of the vector is the magnitude of the frequency and the angle of the vector is the phase. This is the frequency-domain representation of the harmonic.
From trigonometrics it is easy to calculate magnitude and phase form those two values. The magnitude can be calculated from Pythagorean theorem: mag=√(s²+c²) where s and c are measured amplitudes of the respective elements. Similarly phase can be calculated as: phase=arcTan(c/s). You may even find useful to use atan2(c,s) function, which preserves the quadrant.
So we are basically multiplying and then integrating the signal with respective sine and cosine of desired frequency. With complex numbers, there is a shortcut to do that thanks to Euler’s formula:
If we do complex math we can multiply by sine and cosine in a single step. Result will be a complex number – real part will be the magnitude of cosine part and imaginary part will be the magnitude of sine part. We can still calculate the magnitude and phase in very same way. If we do this for every harmonic in the spectrum and plot the complex outputs into graph (3D graph – x = frequency y=real z=imaginary) we have in fact preformed Fourier Transform and we now have the frequency-domain representation of our signal. In the beginning we said that we are analyzing the single period of a harmonic oscillators. Now we know how to measure magnitudes and phases of all the harmonics of the oscillator.
Discrete Fourier Transform formula:
k is the index of the harmonic, zero being DC offset, 1 being the base frequency, 2 being the first overtone etc. N is the length of the analyzed section and x(n) is the n-th element of that section.
Cool fact about Fourier Transform is, that it’s completely reversible. You simply change the sign in the exponent to “+” instead of “-“. This is used for example in additive oscillator, where you define the magnitudes and phases of the harmonics and by inverse Fourier transform you calculate the wave-table for the oscillator. Many times, calculations that are extremely complicated in time-domain are very simple in frequency domain (convolution,decimation and resampling, to name just a few).
Also note, that Fourier transform correlates the input even with frequencies above Nyquist frequency. This makes the testing sine to “alias” and give reflection of the spectrum below Nyquist frequency. At first glance, this doesn’t make any sense, but remember – Fourier Transform is in complex numbers – usually the input is real valued (giving that reflection described above and shown below) however, input may very well be complex-valued giving completely different spectra above Nyquist (yeah, math is crazy).
Off course, the classic (naive) algorithm is a very inefficient way to compute Fourier Transform. Luckily algorithm was invented that can do this much faster, called Fast Fourier Transform. Instead of calculating the frequency domain representations frequency by frequency (which takes enormous calculations, basically N*N complex multiplications and additions per one window) it splits the Time domain signals into single sample bins (Frequency-domain representation of a single value is the single value itself, because it has only the DC component) and then folds them together in frequency domain. That takes only log2(N)*N per window. Flowstone has multiple green primitives that perform FFT and are ready to go. Unfortunately FFT on streams have proven to be more complicated to implement, since nor code nor assembler components have ability to pass arrays as variables. Workarounds are to stream data value by value, or use mem to pointer primitive to pass array pointer to multiple assembler blocks for calculation. Implementation of stream FFT can be seen here by Trogluddite, MyCo and MartinVicanek.
Downside of Fourier Transform is, that it assumes that analyzed signal is periodic and in fact the analyzed part is exactly one period of the signal. That is however rarely true. Often we use Fourier Transform to analyze non-periodic signals, and even with periodic signals we rarely have the opportunity to fit one period exactly into one analysis window. In these cases Fourier Transform may give us slightly unexpected results as shown in the image below. Even when analyzing a simple sine wave, instead of creating one-value peak exactly for that frequency, it shows bell-shaped peak.
This happens almost always, when you use Fourier Transform in analysis of audio. There is a way to reduce these unpleasant artifacts called windowing. Basically you multiply each DFT window with a windowing function that smooths the hard edges on the ends of the window. Many windowing functions were invented each giving slightly different results. Off course, this just makes the frequency spectrum easier to look at and possibly further analyze – it doesn’t really fixes the problem. The real harmonics of the signal lie “in between” the frequency bins (individual harmonics of the analyzed window), so they will always show as bell-shaped peaks.
Analyzing filter Impulse responses vs. analyzing oscillators
As I mentioned before, discrete Fourier transform is based on assumption that the analyzed part of the signal is exactly one period of that signal. That is only half-true though. DFT is also usable, when the analyzed part is the full signal (and only silence is before and after the analyzed section). This is used in analysis of filter frequency and phase response = Fourier transform of the impulse response.
Analyzing filters and oscillators differs in one important thing. When you analyze oscillator and you double the window, frequencies get interleaved with zeros. This happens, because now with double window only even frequencies are present and even frequencies are not – this is obvious, because now period of the oscillator is half the window size.
With filters something completely different happens. Let’s assume that full impulse response of the filter fits into one window (For IIR filters this is theoretically impossible, but in reality, the impulse response decays to zero, so after some point it doesn’t really matter). Now when you double the window size, you just add zeros (or in case of IIR, very small values that don’t really matter) after the IR. When you now preform DFT, the spectrum gets interpolated (becomes twice as dense).
This is interesting fact. It actually means, you can use inverse Fourier transform to generate coefficients for FIR filters. You just pick precision of the response (by picking window size, and by that also the length of generated IR), draw/generate custom frequency and phase response and generate impulse response (which is identical to coefficients) for FIR filter.
OK, hopefully now you are a little bit more familiar with Fourier transform, frequency-domain and what it all means. Have fun incorporating it in your schematics – remember you are limited only by your own fantasy (and power of your computer )!