Fourier Transform in Digital Signal Processing

• Download DiscreteFourierTransform.zip - 126.5 KB
• Download DiscreteFourierTransform.dll.zip - 4.5 KB

Introduction

The frequency analysis is the one of the most popular methods in signal processing. It is a tool for signal decomposition for further filtration, which is in fact separation of signal components from each other.

Although, the process of crossing the border between these two worlds (time and frequency domain) may be considered as advanced case, it is worth doing it, as this other world gives us new opportunities and simplifies many issues.

The article presents implementation of the various versions of calculating Discrete Fourier Transform, starting with definition of Fourier Transform, by reduced calculation algorithm and finishing with Cooley-Tukey method of Fast Fourier Transform.

In addition, it presents the importance of the simplest operations performed on the signal spectrum and their impact on the time domain.

Signals

Signal can be defined as a variability of any physical value, that can be described as a function of a single or multiple arguments. In this article we will be interested in one-dimensional time functions.

In real world, time functions that can be met are placed in continuous domain. However, the development of computer science, caused that analog signal processing became rare. It is much more cost-effective to create, implement and test signal processing algorithms in digital world, then to project and develop analog (electronic) devices.

 

From continuous to discrete domain

In order to receive digital representation of analog signal it needs to be turned into discrete-time domain and quantized.

The Nyquist-Shannon sampling theorem is the link between continuous-time signals and discrete-time signals. This theorem can be also known as The Whittaker-Nyquist-Kotielnikov-Shannon sampling theorem – the choice of the authors names depending on the country in which we talk about this issue. To be above this, we will call it simply the sampling theorem.

The sampling theorem answers the question of how to sample a continuous-time signal to obtain a discrete-time signal, from which you can restore original (continuous-time) signal. According to this statement to obtain a properly sampled discrete-time signal, the sampling frequency must be at least twice of highest frequency that can be observed in original signal.

 

Signal decomposition into frequency domain

It all started in 1807 when the French mathematician and physicist, Joseph Fourier, introduced the trigonometric series decomposition (nowadays know as Fourier series method) to solve the partial differential heat equation in the metal plate. Fourier's idea was to decomposed complicated periodic function into to sum of the simplest oscillating functions - sines ans cosines.

Joseph Fourier

 

If the function f(x) is periodic with period of T and intergrable (its integral is finite) on an interval [x₀, x₀ + T], then it can be transformed into series:

 

where

 

Discrete Fourier Transform

Definition

Fourier series can be named a progenitor of Fourier Transform, which, in case of digital signals (Discrete Fourier Transform), is described with formula:

 

Fourier transformation is reversible and we can return to time domain by calculation:

 

 

In some notations we can observe that divide by N is transterred to inverse calculation - it does not disrupt the calculation unless we apply divide by N both in forward and inverse Fourier Transform.

 

Forward and inverse Discrete Fourier Transform implementations are shown below.

/// <summary>
/// Calculates Discrete Fourier Transform of given digital signal x
/// </summary>
/// <param name="x">Signal x samples values</param>
/// <returns>Fourier Transform of signal x</returns>
public Complex[] DFT(Double[] x)
{
    int N = x.Length; // Number of samples
    Complex[] X = new Complex[N];

    for (int k = 0; k < N; k++)
    {
        X[k] = 0;

        for (int n = 0; n < N; n++)
        {
            X[k] += x[n] * Complex.Exp(-Complex.ImaginaryOne * 2 * Math.PI * (k * n) / Convert.ToDouble(N));
        }

        X[k] = X[k] / N;

    }

    return X;
}



/// <summary>
/// Calculates inverse Discrete Fourier Transform of given spectrum X
/// </summary>
/// <param name="X">Spectrum complex values</param>
/// <returns>Signal samples in time domain</returns>
public Double[] iDFT(Complex[] X)
{
    int N = X.Length; // Number of spectrum elements
    Double[] x = new Double[N];

    for (int n = 0; n < N; n++)
    {
        Complex sum = 0;

        for (int k = 0; k < N; k++)
        {
            sum += X[k] * Complex.Exp(Complex.ImaginaryOne * 2 * Math.PI * (k * n) / Convert.ToDouble(N));
        }

        x[n] = sum.Real; // As a result we expect only real values (if our calculations are correct imaginary values should be equal or close to zero)
    }

    return x;
}

 

The mistake often committed in code implementation of inverse Fourier Transform is to transfer complex value to real by using its magnitude value - in such a case we will receive |x(n)| instead of x(n).

 

Practical example of Discrete Fourier Transform calculated by definition

To verify algorithm let us create a signal that is a sum of two sine waves:

Denote by x₁(n) an original signal and by x₂(n) a disruption signal (noise). The sum of those signals x(n) = x₁(n)+x₂(n) is in that case a disrupted signal.

 

 

   

Fig.1. Original and disruption signals

 

Fig.2. The sum of signals (disrupted signal)

 

As we created our signal from the sum of two sine waves, then according to the Fourier theorem we should receive its frequency image concentrated around two frequencies f₁ and f₂ and also its opposites -f₁ and -f₂.

Fig.3. Amplitude spectrum of disrupted singal.

 

If Fourier theorem is correct, then by removing from spectrum stripes that came from disruption signal, we should receive original signal.

 

 

Fig.4. Spectrum with removed disruption frequency and its inverse DFT result

 

As we can observe in figure 4, our output signal is close to original one and the distortions are caused by additional stripes have arisen as a result of numerical operations. Let's remove them all, and leave only two main stripes.

 

 

Fig.5. Spectrum with only main frequency and its inverse DFT result.

 

The output signal in figure 5 looks like perfect sine wave. But to show minor differences between output and original signal see comparison below.

 

Fig.6. Original and output signal comparison.

 

In figure 3 we cannot see any negative values, so why did we said, that calculated spectrum is concentrated around two frequencies f₁ and f₂ and also its opposites -f₁ and -f₂?

It is because we get use to analyze signals which are defined in range symmetric to the origin. If we denote our disrupted signal x(n) as:

Fig.7. Shifted signal

 

and implement minor change in our algorithm:

/// <summary>
/// Calculates forward Discrete Fourier Transform of given digital signal x with respect to sample numbers
/// </summary>
/// <param name="x">Signal x samples values</param>
/// <param name="SamplesNumbers">Samples numbers vector</param>
/// <returns>Fourier Transform of signal x</returns>
public Complex[] DFT(Double[] x, int[] SamplesNumbers)
{
    int N = x.Length; // Number of samples
    Complex[] X = new Complex[N];

    for (int k = 0; k < N; k++)
    {
        X[k] = 0;

        Double s1 = SamplesNumbers[k]; // Get sample index

        for (int n = 0; n < N; n++)
        {
            Double s2 = SamplesNumbers[n]; // Get sample index
            X[k] += x[n] * Complex.Exp(-Complex.ImaginaryOne * 2 * Math.PI * (s1 * s2) / Convert.ToDouble(N));
        }

        X[k] = X[k] / N;
    }
    return X;
}

 

we will receive amplitude spectrum shown below.

 

Fig.8. Amplitude spectrum of shifted signal

 

The same change needs to be implemented in inverse DFT method.

/// <summary>
/// Calculates inverse Discrete Fourier Transform of given spectrum X with respect to sample numbers
/// </summary>
/// <param name="X">Spectrum complex values</param>
/// <param name="SamplesNumbers">Samples numbers vector</param>
/// <returns>Signal samples in time domain</returns>
public Double[] iDFT(Complex[] X, int[] SamplesNumbers)
{
    int N = X.Length; // Number of spectrum elements
    Double[] x = new Double[N];

    for (int n = 0; n < N; n++)
    {
        Complex sum = 0;
        Double s1 = SamplesNumbers[n]; // Get sample index

        for (int k = 0; k < N; k++)
        {
            Double s2 = SamplesNumbers[k];  // Get sample index
            sum += X[k] * Complex.Exp(Complex.ImaginaryOne * 2 * Math.PI * (s1 * s2) / Convert.ToDouble(N));
        }

        x[n] = sum.Real; // As a result we expect only real values (if our calculations are correct imaginary values should be equal or close to zero)
    }

    return x;
}

 

Notice, that by shifting signal to symmetrical position to the orgin, we calculated shifted spectrum:

 

 

The same effect can be achieved, if we change original signal using formula

 

 

In this case, there is no need to store the actual numbers of samples. For this use the method is shown below.

 

/// <summary>
/// Shifts spectrum by changing original signal before Fourier Transform
/// </summary>
/// <param name="x">Original signal</param>
/// <returns>Shifted signal</returns>
public Double[] Shift(Double[] x)
{
    int N = x.Length;

    for (int i = 0; i < N; i++)
    {
        x[i] = (int)Math.Pow(-1, i) * x[i];
    }

    return x;
}

 

Discrete Fourier Transform with a reduced number of calculations

Let us look again at the forward Discrete Fourier Transform definition:

 

 

Denote by W_N auxiliary variable:

 

 

Then our formula takes the form

 

 

Also we can notice that if k = 0 or n = 0 we receive:

 

 

Accordingly, the designation of DTF can be written in matrix form:

 

 

It follows that, to determine the Dicrete Fourier Transform of N-samples signal, we have to perform N² multiplication operations and determine the (N-1)² coefficients in matrix W_N. Remind that multiplication opertion is one of the most time-consuming opertion for processors.

Let us return to factor W_N and consider N=4 samples signal x(n). To calculate X(k) we need to find values of matrix coefficients:

We have:

where

is the Euler's formula.

 

Next values are

 

Notice that

 

Similar situation can be noticed for N=8:

 

 

that provides to:

 

In general, for N-samples signal we can write:

 

This observation leads us to the conclusion that in determining the DFT we calculate only N/2 - 1 spectral coefficients, and another we can get from the previous ones - it is a property of symmetry of the spectrum. Due to this observation four times reduced the number of operations within the algorithm. This observation is true only if we deal with the physical signals, that is signals with real values.

/// <summary>
/// Calculates forward Discrete Fourier Transform of given digital signal x with reduced number of multiplications
/// </summary>
/// <param name="x">Signal x samples values</param>
/// <returns>Fourier Transform of signal x</returns>
public Complex[] DFT2(Double[] x)
{
    int N = x.Length; // Number of samples
    Complex[] X = new Complex[N];

    // Zero-index element calculation
    X[0] = x.Sum() / N;

    for (int k = 1; k < N/2 + 1; k++)
    {
        X[k] = 0;

        for (int n = 0; n < N; n++)
        {
            X[k] += x[n] * Complex.Exp(-Complex.ImaginaryOne * 2 * Math.PI * (k * n) / Convert.ToDouble(N));
        }

        X[k] = X[k] / N;

        // Miss additional calculation for center element
        if (k != N / 2)
        {
            X[N - k] = new Complex(X[k].Real, -X[k].Imaginary);
        }
    }

    return X;
}

/// <summary>
/// Calculates forward Discrete Fourier Transform of given digital signal x with respect to sample numbers and reduced number of multiplications
/// </summary>
/// <param name="x">Signal x samples values</param>
/// <param name="SamplesNumbers">Samples numbers vector</param>
/// <returns>Fourier Transform of signal x</returns>
public Complex[] DFT2(Double[] x, int[] SamplesNumbers)
{
    int N = x.Length; // Number of samples
    Complex[] X = new Complex[N];

    for (int k = 0; k < N / 2 + 1; k++)
    {
        X[k] = 0;
        Double s1 = SamplesNumbers[k];

        for (int n = 0; n < N; n++)
        {
            Double s2 = SamplesNumbers[n];
            X[k] += x[n] * Complex.Exp(-Complex.ImaginaryOne * 2 * Math.PI * (s1 * s2) / Convert.ToDouble(N));
        }

        X[k] = X[k] / N;

        // Miss additional calculation for center element
        if (k != N/2)
        {
            X[N - k - 1] = new Complex(X[k].Real, -X[k].Imaginary);
        }
    }

    return X;
}

 

Examples:

 

Fig.9. Amplitude spectrum of signal from fig.2. designated with DFT with reduced number of calculation

 

Fig.10. Amplitude spectrum of signal from fig.7. designated with DFT with reduced number of calculation

 

Fast Fourier Transform

Fourier Transform formula can be splitted into:

where

and

therefore:

that can be written in a simpler form, that separates the calculations performed for the even samples, from those performed for the odd ones:

Accordingly, referring to our previous calculations for the case N = 8 can be written as:

 

which can be presented as graph:

 

Fig.11. DFT decomposition into even and odd samples calculation.

 

Note that the elements: 

   and   

 

can be calculated in the same way: elements x(0), x(4), x(1) and x(5) are in their diagrams the even indices, whereas x(2), x(6), x(3) and x(7) are the odd ones. Combining them in the ordered pairs, the above graphs can be redraw:

Finally, in this example:

Fig.12. Fast Fourier Transform (RADIX-2) for N = 8 signal.

 

Designating DFT definition for N = 8 signal we had to perform 64 operations of multiplication, but due to the above observation, we made them only 12. The greater the number of samples in the signal, the greater benefit in calculations number with the use of the FFT algorithm, since the amount of calculations for N-smaple signal can be determined as:

The basic element of such a calculation is called "butterfly":

 

where

 

Successive decomposition stages for smaller, down to the basic "butterflies" may be carried out only on the assumption that N is a power of 2.

 

Fast Fourier Transform in practice

 

Proper number of samples

The easiest way to processed signal have the number of samples which is a power of 2, simply fill in the missing samples of zeros. Depending on the signal structure, additions can be made as right-hand, left-hand or both sides.

 

Fig.13. Complemented signal from fig.2

 

Sorting samples

To connect consecutive signal samples in the respective pairs, we are using preliminary sorting algorithms. However, the FFT algorithm can be performed on samples unsorted, and the sorting operation (in the same way) can be performed on the samples of the spectrum.

In the hardware spectrum analyzers and low-level programming languages sorting samples it is much easier, because the operation is based on a simple inversion of bits in the binary representation of the index sample.

In our N = 8 signal we can notice that:

Method implementation in general case:

/// <summary>
/// Input data sorting using bit reversal method
/// </summary>
/// <param name="x">Original input signal samples</param>
/// <returns>Sorted signal samples</returns>
public Double[] FFTDataSort(Double[] x)
{
    int N = x.Length; // signal length

    if (Math.Log(N, 2) % 1 != 0)
    {
        throw new Exception("Number of samples in signal must be power of 2");
    }

    Double[] y = new Double[N]; // output (sorted) vector

    int BitsCount = (int)Math.Log(N, 2); // maximum number of bits in index binary representation
    
    for (int n = 0; n < N; n++)
    {
        string bin = Convert.ToString(n, 2).PadLeft(BitsCount, '0'); // index binary representation
        StringBuilder reflection = new StringBuilder(new string('0', bin.Length));

        for (int i = 0; i < bin.Length; i++)
        {
            reflection[bin.Length - i - 1] = bin[i]; // binary reflection
        }

        int number = Convert.ToInt32(reflection.ToString(),2); // new index

        y[number] = x[n];

    }

    return y;

}

 

This operation can, however, be performed much faster, without need to convert the index number to the binary form:

/// <summary>
/// Input data sorting using index shift method
/// </summary>
/// <param name="x">Original input signal samples</param>
/// <returns>Sorted signal samples</returns>
public Double[] FFTDataSort2(Double[] x)
{
    int N = x.Length; // signal length

    if (Math.Log(N, 2) % 1 != 0)
    {
        throw new Exception("Number of samples in signal must be power of 2");
    }

    int m = 1; // original index
    
    for (int n = 0; n < N - 1; n++) // n - new index
    {
        if (n < m)
        {
            Double T = x[m-1];
            x[m-1] = x[n];
            x[n] = T;
        }

        int s = N / 2; // shift operator

        while (s < m)
        {
            m = m - s;
            s = s / 2;
        }

        m = m + s;
    }

    return x;

}

 

Fig.14. Sorted signal samples (signal from fig.2. complemented at the right side)

 

Fast Fourier Transfom algorithm

FFT algorithm implementation is shown below.

/// <summary>
/// Calculates forward Fast Fourier Transform of given digital signal x
/// </summary>
/// <param name="x">Signal x samples values</param>
/// <returns>Fourier Transform of signal x</returns>
public Complex[] FFT(Double[] x)
{
    int N = x.Length; // Number of samples

    if (Math.Log(N, 2) % 1 != 0)
    {
        throw new Exception("Number of samples in signal must be power of 2");
    }


    Complex[] X = new Complex[N]; // Signal specturm

    // Rewrite real signal values to calculated spectrum
    for (int i = 0; i < N; i++)
    {
        X[i] = new Complex(x[i], 0);
    }

    int S = (int)Math.Log(N, 2); // Number of calculation stages

    for (int s = 1; s < S + 1; s++) // s - stage number
    {
        int BW = (int)Math.Pow(2, s - 1); // the width of butterfly
        int Blocks = (int)(Convert.ToDouble(N) / Math.Pow(2, s)); // Number of blocks in stage
        int BFlyBlocks = (int)Math.Pow(2, s - 1); // Number of butterflies in block
        int BlocksDist = (int)Math.Pow(2, s); // Distnace between blocks

        Complex W = Complex.Exp(-Complex.ImaginaryOne * 2 * Math.PI / Math.Pow(2, s)); // Fourier Transform kernel

        for (int b = 1; b < Blocks + 1; b++)
        {
            for (int m = 1; m < BFlyBlocks + 1; m++)
            {
                int itop = (b - 1) * BlocksDist + m; // top butterfly index
                int ibottom = itop + BW; // bottom butterfly index

                Complex Xtop = X[itop-1]; // top element -> X(k)
                Complex Xbottom = X[ibottom-1] * Complex.Pow(W, m - 1); // bottom element -> X(k + N/2)

                // Butterfly final calculation
                X[itop-1] = Xtop + Xbottom;
                X[ibottom-1] = Xtop - Xbottom;
            }
        }
    }

    // Spectrum scaling
    for (int i = 0; i < N; i++)
    {
        X[i] = X[i] / Convert.ToDouble(N);
    }

    return X;
}

 

Fig.15. FFT - Amplitude spectrum of signal with sorted samples (from fig.14)

 

Fig.16. FFT - Amplitude spectrum of shifted signal with sorted samples (from fig.14)

 

Inverse Fast Fourier Transform can be implemented by minor change in forward FFT method:

/// <summary>
/// Calculates inverse Fast Fourier Transform of given spectrum
/// </summary>
/// <param name="X">Spectrum values</param>
/// <returns>Signal samples</returns>
public Double[] iFFT(Complex[] X)
{
    int N = X.Length; // Number of samples
    Double[] x = new Double[N];
    int E = (int)Math.Log(N, 2);

    for (int e = 1; e < E + 1; e++)
    {
        int SM = (int)Math.Pow(2, e - 1);
        int LB = (int)(Convert.ToDouble(N) / Math.Pow(2, e));
        int LMB = (int)Math.Pow(2, e - 1);
        int OMB = (int)Math.Pow(2, e);

        Complex W = Complex.Exp(Complex.ImaginaryOne * 2 * Math.PI / Math.Pow(2, e)); // changed sign - our minor change

        for (int b = 1; b < LB + 1; b++)
        {
            for (int m = 1; m < LMB + 1; m++)
            {
                int g = (b - 1) * OMB + m;
                int d = g + SM;

                Complex xgora = X[g - 1];
                Complex xdol = X[d - 1] * Complex.Pow(W, m - 1);

                X[g - 1] = xgora + xdol;
                X[d - 1] = xgora - xdol;
            }
        }
    }

    for (int i = 0; i < N; i++)
    {
        x[i] = X[i].Real;
    }

    return x;
}

 

Fig.17. Signal restored from spectrum from fig.15 using inverse Fast Fourier Transform.

 

As it can be noticed, the signal complemented with additional zero samples is not without effect on the spectrum - at figures 15 and 16 additional stripes can be noticed. Nevertheless, while we keep them for futher calculations, we are able to restore the original signal.

However, let us see if we can still separate the signal from noise and restore it using FFT methods.

 

Firstly, remove stripes that are suspected of being a disturbance:

 

Fig.18. Filtred signal spectrum

 

Secondly, calculate inverse FFT

 

Fig.19. Inverse FFT - restored signal

 

Thirdly, cut additional samples

 

Fig.20. Signal without additional samples.

 

At the end, compare original (undisturbed) signal with our filtration result.

 

Fig.21. Original and output signal comparison (Fast Fourier Transform).

 

Using the code

All signal processing algorithms presented in this article are implemented in attached source code. To start working with DiscreteFourierTransform library simply declare its namespace:

using DiscreteFourierTransform;

create new object:

DiscreteFourierTransform.DiscreteFourierTransform dft = new DiscreteFourierTransform.DiscreteFourierTransform();

and start your adventure with digital signal processing in frequency domain.

 

Library methods and applications

Test signal

Disturbed sine signal and samples indexes are included in dll:

Double[] signal = dft.testsignal();
int[] n = dft.testsignalindexes();

 

Discrete Fourier Transform

// --- Discrete Fourier Transform by definition

Complex[] X = dft.DFT(signal); // discrete fourier transform by definition

Double[] A = dft.Amplitude(X); // amplitude spectrum from X

Complex[] Xa = dft.DFT(signal, n); //discrete fourier transform by definition with samples numbers

Double[] signal2 = dft.Shift(signal); // Signal shift method

Complex[] Xb = dft.DFT(signal2); // discrete fourier transform by definition of shifted signal

// Compare Xa and Xb and see if we were correct!

 

Inverse Discrete Fourier Transform

// --- Inverse Discrete Fourier Transform by definition

Double[] x = dft.iDFT(X); // inverse discrete fourier transform by definition

Double[] xa = dft.iDFT(Xa, n); // inverse discrete fourier transform by definition (with samples numbers)

Double[] xb = dft.iDFT(Xb); // inverse discrete fourier transform by definition of shifted signal

// Compare x, xa and xb. Any other operation needs to be performed at xb to restore original?

 

Discrete Fourier Transform with reduced number of calculation

// --- Discrete Fourier Transform with reduced number of calculation

Complex[] X2 = dft.DFT2(signal); // discrete fourier transform with reduced number of calculation

Complex[] X2a = dft.DFT2(signal, n); // discrete fourier transform with reduced number of calculation with samples numbers

Complex[] X2b = dft.DFT2(dft.Shift(signal), n); // discrete fourier transform with reduced number of calculation (shifted signal)

 

Fast Fourier Transform

Signal pack with zeros

// Signal pack with zeros

Double[] signalpack = dft.SignalPackRight(signal); // Complement at right to N = nearest power of 2

Double[] signalpack2 = dft.SignalPackLeft(signal); // Complement at left to N = nearest power of 2

Double[] signalpack3 = dft.SignalPackBoth(signal); // Complement at both sides to N = nearest power of 2

 

Samples sorting

// Samples sorting

Double[] signal3 = dft.FFTDataSort(signalpack); // using bit reverse

Double[] signal4 = dft.FFTDataSort2(signalpack); // using shift method

// Compare signal3 and signal4 - should be equal!

 

Fast Fourier Transform

Complex[] X3 = dft.FFT(signal3); // fast fourier transform

Complex[] X3a = dft.FFT(dft.Shift(signal3)); // shifted spectrum - fast fourier transform

 

Inverse Fast Fourier Transform

Spectrum stripes sorting

Complex[] X3s = dft.FFTDataSort(X3); // spectrum sort using bit reverse

Complex[] X3as = dft.FFTDataSort2(X3); // spectrum sort using shift method

 

Inverse Fast Fourier Transform

Double[] x3 = dft.iFFT(X3s); // inverse fast fourier transform

// ----------------- Compare x3 with signal!!! ---------------