Discrete Fourier Transform

Introduction:

Viewing the frequency spectra of signals is useful for filtering of signals and for signal extraction. To convert a signal to the frequency domain, numerous transforms and algorithms exist. Frequency spectrum conversions for signals include the Fourier Transform(FT), Discrete Fourier Transform(DFT), Fast Fourier Transform(FFT), Laplace Transform, Z-Transform…ect. With a wide variety of conversions, it is easy to wonder which is the best and most useful transformation? Naturally, the answer is that each transformation has its own advantages and disadvantages. This article will present how to use the Discrete Fourier Transform to analyze data. The advantages of the DFT is that it is very easy to program and consequently is easy to implement in any coding language. The disadvantage is the computational time required but this is becoming less and less of an issue in the presence of modern computing. An equation for a DFT is presented below:

So how does it work? In the continuous time domain, integrating and multiplying a signal by a complex exponential ($e^{-j 2 \pi f \cdot t}$) decomposes a signal into its frequency components. In discrete time, integration becomes a summation. In the above equation, summing variable $k$ integrates through various frequencies while summing variable $n$ iterates through the values of the discrete signal $x[n]$.

Basic Example:

To demonstrate this process, lets look at a simple test signal $x[n] = cos(10 \pi t)$ and plot it in matlab/octave.

clc;clear;
Step = 0.001;
t= 0:Step:1;
x = cos(10*pi*t);
plot(t,x,'linewidth',2)
axis([0 1 -1.2 1.2])
title('cos(10\pit)')
grid on

Next we will perform the DFT on the signal in order to observe its frequency components. Naturally a cosine is of a pure singular frequency component. In this case, the cosine wave $x[n]=cos(10 \pi n)$ has a frequency of $5 Hz$. Note that frequency of $cos(k \cdot x)$ can be calculated with $f = k / 2 \pi = 5$. The following image is the DFT waveform:

Note that in this waveform, the DFT is not simply $\delta ( k – 5)$ as we would expect. Reasons for this includes that the DFT is an approximation since cosine is not of infinite resolution, also that the DFT is performed only using $k$ amount of spectra, and that the signal sampled is not of infinite period. Consequently, there are artifacts of many other signals than just 5 Hz. None the less, this approximation clearly shows that 5 Hz is the dominate frequency spectra of the signal. In regards to coding the DFT, built in matlab / octave functions can be used. However, to gain a better understanding of the algorithm, I implemented it in matlab / octave with the following code:

function [DFT,DFTX] = DFTP(X,TimeStep,Points,maxFreq)
N=length(X);
j = (-1)^(1/2);

if nargin < 4 % if max freq is passed in, use it, else use alis freq
    AliasFreq = 1/TimeStep/2;
else
    AliasFreq = maxFreq;
end

Total_Time = TimeStep*length(X);
DFTX = 1:Points;
DFTX = DFTX*AliasFreq/Points;

DFT = zeros(Points,1);

for n=1:N
    for k=1:Points
        DFT(k)=DFT(k)+ X(n)*exp(-j*2*pi*k*AliasFreq*Total_Time/(N*Points)*n);
    end
end

DFT = abs(DFT); % Remove imaginary component, look at magnitude 
DFT = DFT./sum(DFT); %represent signal magnitude as PDF

end

This function allows for the user to specify the resolution of the DFT with $Points$. It also allows a user to specify the time step of the signal as well as the max frequency. If a max frequency is not specified, the anti alias frequency of $f = 1 / (TimeStep * 2)$ is used. This function can be called and ploted with the following code:

clc;clear;
Step = 0.001;
t= 0:Step:3;
x = cos(10*pi*t);
[DFTY,DFTX] = DFTP(x,Step,1000,200);
plot(DFTX,DFTY);

More Complex Example:

Next, to demonstrate to help demonstrate the usefulness of the DFT, a more complicated signal can be plotted and white noise $e(t)$ is added:

$$x[n]=sin\big( 150(2\pi\cdot n)\big) + sin\big( 40(2\pi\cdot n)\big) + sin\big( 5(2\pi\cdot n)\big) + 4\cdot e(n) $$

When the above signal is plotted and the DFT is performed, we observe the following:

Clearly the dominate spectra are at $f=5,40,150Hz$ with roughly even distribution.

Making sense of DFT Magnitude:

Note that my DFT function’s height is in terms of a probability mass function(pdf, $f_n[n]$) for discrete time. This is similar to the probability density function for continuous time. The probability mass function is defined by the following characteristic.

$$\sum^{N}_{k=0}f_k[k]=1 \\ P(k=X) = f_k[X] $$

Note the height of a DFT is somewhat arbitrary which is why I created a pdf. A one-normalization would have worked as well but this would have been poor at comparing signals. A demonstration of comparing DFT pdfs can be seen below:

Here the more pure signal $y[n]$ is of higher magnitude than the combined signal $x[n]$. This difference is useful when calculating the cross correlation or cross co-variance. This difference would have been nullified had the signals been normalized. For reference, the complete source code of the above plot is included below:

clc;clear;
figure(1)
subplot(2,1,1)
Step = 0.001;
t= 0:Step:3;
xx = cos(200*pi*t)+1*rand(1,length(t))-0.5;
x = 1*sin(150*(2*pi*t)) + 1*sin(40*(2*pi*t))+sin(5*(2*pi*t))+1*rand(1,length(t))-0.5;
plot(t,x,t,xx,'linewidth',1)
axis([0 0.4 -3 3.5])
title('cos(10\pit)')
legend('x[n]','y[n]')
grid on
title('x[n]=sin(150(2\pi\cdotn)) + sin(40(2\pi\cdotn)) + sin(5(2\pi\cdotn)) + \e(n)\newliney[n]=cos(100(2\pi\cdotn))+e(n)')
%figure(2)
subplot(2,1,2)
[DFTY,DFTX] = DFTP(x,Step,1000,200);
[DFTYY,DFTXX] = DFTP(xx,Step,1000,200);
plot(DFTX,DFTY,DFTXX,DFTYY,'linewidth',1);
title('DFT(x[n]), DFT(y[n])')
xlabel('Hz')
ylabel('pdf')
ylim([0,0.08])
legend('DFT(x[n])','DFT(y[n])')
grid on

References:

https://en.wikipedia.org/wiki/Discrete_Fourier_transform
Many thanks to Dr. Josse of Marquette University for his fundamental probability course.