Tag Archives: DFT

Fast Fourier Transform – Code

Fundamentals, , , ,

Most of us have used the FFT routine in MATLAB. This routine has become increasingly important in simulation of communication systems as it is being used in Orthogonal Frequency Division Multiplexing (OFDM) which is employed in 4G technologies like LTE and WiMAX. We would not go into the theoretical details of the FFT, rather, we would produce the MATLAB code for it and leave the theoretical discussion for a later time.

The underlying technique of the FFT algorithm is to divide a big problem into several smaller problems which are much easier to solve and then combine the results in the end.

%%%%% INITIALIZATION %%%%%
clear all
close all

fm=100;
fs=1000;
Ts=1/fs;
nn=1024;
isign=1;

t=0:Ts:(nn-1)*Ts;

x_c=exp(1i*2*pi*fm*t);
x(1:2:2*nn-1)=imag(x_c);
x(2:2:2*nn)=real(x_c);

%%%%% BIT REVERSAL %%%%%
n=2*nn;
j=1;
for i=1:2:n-1
    if(j>i)
        temp=x(j);
        x(j)=x(i);
        x(i)=temp;

        temp=x(j+1);
        x(j+1)=x(i+1);
        x(i+1)=temp;
    end
    m=nn;
    while (m >= 2 && j > m)
        j=j-m;
        m=m/2;
    end
    j=j+m;
end

%%%%% DANIELSON-LANCZOS ALGO %%%%%
mmax=2;
while (n > mmax)
    istep=2*mmax;
    theta=isign*(2*pi/mmax);
    wtemp=sin(0.5*theta);
    wpr=-2*wtemp*wtemp;
    wpi=sin(theta);
    wr=1.0;
    wi=0.0;
    for m=1:2:mmax-1
        for i=m:istep:n
            j=i+mmax;
            tempr=wr*x(j)-wi*x(j+1);
            tempi=wr*x(j+1)+wi*x(j);
            x(j)=x(i)-tempr;
            x(j+1)=x(i+1)-tempi;
            x(i)=x(i)+tempr;
            x(i+1)=x(i+1)+tempi;
        end
        wtemp=wr;
        wr=wr*wpr-wi*wpi+wr;
        wi=wi*wpr+wtemp*wpi+wi;
    end
mmax=istep;
end
F=x(2:2:n)+1i*x(1:2:n-1);
plot((0:fs/(nn-1):fs),abs(F))

In the above example we have calculated an ‘nn’ point complex FFT of an ‘nn’ point complex time domain signal. The algorithm takes the input in a special arrangement where the ‘nn’ point complex input signal is converted into ‘2*nn’ real sequence where the imaginary components are placed in odd elements and real components are placed in even elements of the input sequence. A similar arrangement works for the output sequence.

Shown below is the FFT of a complex exponential with a frequency of 100 Hz. The plot is shown from 0Hz to 1000 Hz which is the sampling frequency. A signal with multiple frequencies would have to be passed through a Low Pass Filter (LPF) so that the signal components above 500 Hz (fs/2) are filtered out. When the FFT of a real signal is performed an image frequency is produced between 500 Hz to 1000 Hz.

Fast Fourier Transform of a Complex Exponential
Fast Fourier Transform of a Complex Exponential

Here we have discussed the case of complex input sequence. Simplifications can be made for a real sequence or for special signals such as pure sine and cosine waves. We will discuss these in later posts.

 [1] Numerical Recipes in C, The Art of Scientific Computing, Cambridge University Press.

Detecting Sinusoids in Noise

Fundamentals, , , ,

We have previously discussed the problem of detecting two closely spaced sinusoids using the Discrete Fourier Transform (DFT). We assumed that the data set we got was pure i.e. there was no noise. However, in reality this is seldom the case. There is always some noise, corrupting the signal. Let us now see how it effects the detection problem.

We consider Additive White Gaussian Noise (AWGN) as the corrupting source. The noise power is set equal to the power of the two sinusoids i.e. we have an SNR of 0 dB. This is quite a severe case, the noise power is usually a few dB below the signal power. We are also bounded by the number of samples, N=64, giving us a resolution of 15.87 Hz.

clear all
close all

fm1=100;
fm2=120;
fs=1000;
Ts=1/fs;
N=64;
t=0:Ts:(N-1)*Ts;

x1=sqrt(2)*cos(2*pi*fm1*t);
x2=sqrt(2)*cos(2*pi*fm2*t);
wn=randn(1,N);
x=x1+x2+wn;
W=exp(-j*2*pi/N);

n=0:N-1;
k=0:N-1;

X=x*(W.^(n'*k));

plot(k/N,20*log10(abs(X)))
xlabel ('Normalized Frequency')
ylabel ('X(dB)')

The data is plotted on a logarithmic scale so that we can compare the signal and noise  power levels.

DFT of Two Tones in Noise
DFT of Two Tones in Noise

It is observed that we still have two peaks around the required frequency bins but there are also a number of false peaks. These peaks are around 10 dB lower than the signal peaks and should not cause a false detection. So the signal to noise ratio of 0 dB in the time domain is translated to a signal to noise ratio of about 10 dB in the frequency domain (this can be realized using an appropriate filter).

Resolution of Discrete Fourier Transform

Fundamentals, , ,

In the previous post we had introduced the Discrete Fourier Transform (DFT) as a method to perform the spectral analysis of a time domain signal. We now discuss an important property of the DFT, its spectral resolution i.e. its ability to resolve two signals with similar spectral content.

Initially one might think that increasing the sampling frequency would increase the spectral resolution but this totally incorrect. In fact if the sampling frequency is increased, keeping the number of time domain samples to be the same, the resolution actually decreases. So how do we calculate the spectral resolution. One simple way is to calculate the difference between two frequency bins as fs/(N-1) or 1/(N-1)Ts. Simply put the resolution in the frequency domain is the inverse of the sample length in the time domain.

So let us now calculate the DFT of two closely spaced sine waves keeping the sampling frequency to be the same and changing the number of time domain samples (only the result for N=64 shown here). We again list down the code used to calculate the DFT.

clear all
close all

fm1=100;
fm2=120;
fs=1000;
Ts=1/fs;
N=64;
t=0:Ts:(N-1)*Ts;

x1=cos(2*pi*fm1*t);
x2=cos(2*pi*fm2*t);
x=x1+x2;
W=exp(-j*2*pi/N);

n=0:N-1;
k=0:N-1;

X=x*(W.^(n'*k));

plot(k/N,abs(X))
xlabel ('Normalized Frequency')
ylabel ('X')

We first do some quick math to find the spectral resolution.

 fs/(N-1)=1000/(64-1)=15.87 Hz

Now the two tones are space 20 Hz apart (100 Hz and 120 Hz), so we can predict that the two tones would be detected successfully. The result of the DFT operation on the composite signal is shown below.

DFT of Two Tones
DFT of Two Tones

It is observed that although the two tones are detected, they are not exactly at the desired frequencies (0.10 and 0.12). Secondly the amplitude of the two tones is different although the time domain signals had equal amplitude. Both these phenomenon are due to the fact that we only have a limited number of frequency bins (N=64) due to which the resulting spectrum is only an estimate of the true spectrum.

There are better techniques than DFT to separate two closely spaced sinusoids and these are known as super resolution spectral techniques and would be discussed some other time.

Discrete Fourier Transform

Fundamentals, ,

Discrete Fourier Transform or DFT is a mathematical operation that transforms a time domain signal to frequency domain. It is usually implemented using the Fast Fourier Transform (FFT). The computational complexity of the DFT is N2 whereas its (N)log2N for the FFT, where N is the number of samples of the the time domain signal. Mathematically, the DFT is written as

and this operation can be easily implemented in MATLAB as shown below.

clear all
close all

fm=100;
fs=1000;
Ts=1/fs;
N=1024;
t=0:Ts:(N-1)*Ts;

x=cos(2*pi*fm*t);
W=exp(-j*2*pi/N);

n=0:N-1;
for k=0:N-1
        X(k+1)=sum(x.*(W.^(k*n)));
end

plot((0:k)/N,abs(X))
xlabel ('Normalized Frequency')
ylabel ('X')

The DFT of the cosine function with frequency of 100Hz is shown below. It must be noted that the frequency in the figure below is normalized by the sampling frequency. So the peaks in the graph occur at 0.1*1000=100Hz and 0.9*1000=900Hz (image frequency). Actually if you zoom in you will see that the peaks are not exactly at 0.1 and 0.9 but  are slightly offset as the frequency bins are not located at exactly those frequencies.

DFT of a cosine wave
DFT of a cosine wave

In case of multiple sinusoids the resolution of the DFT becomes important. Higher the sample size (duration of the time domain signal) higher is the resolution of the DFT. As stated earlier the DFT is a computationally complex operation and usually the Fast Fourier Transform (FFT) is used to compute the frequency domain behavior. We will discuss this in the following posts.

A simple way to expedite the process of DFT calculation is to use matrix manipulation instead of a “for loop”.This is shown below.

clear all
close all

fm=100;
fs=1000;
Ts=1/fs;
N=128;
t=0:Ts:(N-1)*Ts;

x=cos(2*pi*fm*t);
W=exp(-j*2*pi/N);

n=0:N-1;
k=0:N-1;

X=x*(W.^(n'*k));

plot(k/N,abs(X))
xlabel ('Normalized Frequency')
ylabel ('X')

Although matrix manipulation makes for a nice and clean code, it was found that there was no improvement in the computation time. Both the “for loop” code and the “matrix multiplication” code took about 0.35 seconds to execute. However, increasing the sample size from 128 to 1024 results in significant better computation time for the latter scheme.