Demystifying the Fourier Transform

The Fourier Transform is often used in Communication and Signal Processing to find the spectral content of a time-domain signal. The most common example is that of a sinusoid in the time domain, resulting in a sharply peaked signal in the frequency domain, also known as a delta function. A rectangular pulse in the time domain has a more complicated frequency domain equivalent, a sinc function. A rectangular pulse may be thought of as a combination of many sinusoids, hence its frequency domain equivalent is not that straightforward.

Continue reading Demystifying the Fourier Transform

A Comparison of FFT, MUSIC and ESPRIT Methods of Frequency Estimation

As discussed in previous posts it is frequently required in communications and signal processing to estimate the frequency of a signal embedded in noise and interference. The problem becomes more complicated when the number of observations (samples) is quite limited. Typically, the resolution in the frequency domain is inversely proportional to the window size in the time domain. Sometimes the signal is composed of multiple sinusoids where the frequency of each needs to be estimated separately. Simple techniques such as Zero Crossing Estimator fail in such a scenario.  Even some advanced techniques such as MATLAB function “pwelch” fail to distinguish closely spaced sinusoids.

Continue reading A Comparison of FFT, MUSIC and ESPRIT Methods of Frequency Estimation

Frequency Estimation Using Zero Crossing Method

A sinusoidal signal is the most fundamental type of signal that exists in communication systems, power systems, navigation systems etc. It is controlled by three parameters which are the amplitude, phase and frequency. The last two, that is phase and frequency, are interconnected. As discussed in my previous post Instantaneous Frequency (IF) is nothing but the rate of change of phase. This can be mathematically described as:

IF=Δφ/Δt

Continue reading Frequency Estimation Using Zero Crossing Method

Fourier Transform in Python

Fast Fourier Transform or FFT is a powerful tool to visualize a signal in the frequency domain. Shown below is the FFT of a signal (press the play button) composed of four sinusoids at frequencies of 50Hz, 100Hz, 200Hz and 400Hz. The sampling frequency is set at 1000Hz, more than twice the maximum frequency of the composite signal. The amplitude of the sinusoids diminishes with increasing frequency and this is also reflected in the frequency domain. Play with the code, decrease the frequencies, increase the power, visualize the FFT output on logarithmic scale!

Fast Fourier Transform – Code

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.

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

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.

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).