An article dedicated to learn and develop some basic approaches related to Discrete Fourier Transform.

First of all, I would encourage tp read a great article there. The great simplified approach and insights were provided by that link.

Let’s go through some key point and gain better understanding and develop some useful skills.

Correlation is a first key concept there, it can be described as the sum of the product of two signals:

\sum_{i=0}^{N}x(i)y(i)

We are talking there about the simplest correlation form, without any cross correlation equations involved.

The correlation can be positive, or negative if the signal is the same but simply out of the phase, or if signals are completely different correlation will have extremely small values.

So the key point in whole discrete transformation just a simple sweep – you take the point a the certain frequency (which you know can be represented by Euler equation) multiply by your signal, check the correlation, then continue up to the end.

The easiest way to learn is to look at some simple examples, let’s say we have sine and cosine at frequency equal to 7Hz:

And next there is a simple Octave piece of the code:

clf
clear
#input data variables
step = 1/256;
tstop = 2;
f1 = 7;

t = 0:step:tstop-step;

for i = 1:length(t)
  rand_coeff1(i) = (rand()-0.5:0.5)/3;
  rand_coeff2(i) = (rand()-0.5:0.5)/3;
end
y = sin(2*pi*f1*t);

ytry=e.^(j*2*pi*7*t);
y_product1 = y.*ytry;

plot(t,y)
hold on
grid on

correlation1 = sum(imag(y_product1));
correlation2 = sum(real(y_product1));


y = sin(2*pi*f1*t+pi/2);
plot(t,y)
ytry=e.^(j*2*pi*7*t);
y_product1 = y.*ytry;

correlation3 = sum(imag(y_product1));
correlation4 = sum(real(y_product1));

Here, in that code, one important thing should pop up – we should not forget to pull the phase information from these products.

It could be clearly seen that correlation1 is indicating that signal which has been probed is the sine, while correlation4 indicates the cosine.

OK, so we can now distinguish sine from cosine, let’s make our example a little more complex – now we have two signals – at 7 and at 8Hz:

clf
clear
#input data variables
step = 1/256;
tstop = 2;
f1 = 7;

t = 0:step:tstop-step;

N = 16;

y = sin(2*pi*f1*t)+sin(2*pi*(f1+1)*t);

#ytry=e.^(j*2*pi*1*t);
#y_product1 = y.*ytry;

#plot(t,y)
##hold on
#grid on


for i=1:1:N
  ytry=e.^(j*2*pi*i*t);
  y_product = y.*ytry;
 
  correlation_real(i) = sum(real(y_product));
  correlation_imag(i) = sum(imag(y_product));
  correlation_abs(i) = correlation_real(i) + correlation_imag(i);
  plot(i,correlation_abs(i),'marker','o');
  hold on;  
  grid on;
endfor

Checking the result:

Two peaks clearly seen at 7 and 8, looking into that picture I had couple questions – Why are these numbers so big? How to define proper N? What is the limits? How to make plot with bars instead of markers 🙂

Here some more realistic function, this time lets do it with numbers, so we will have some kind of an anсhor for us:

y = 2 + 5cos(2 \pi2 t)+3sin(2\pi3t)+2cos(2\pi4t)

And the sampled one curve is a set of values sampled with Fs = 8Hz, here I took 8 samples, which means that resolution will be equal to 1Hz.

To answer question why numbers are that big, it could be seen that it just corresponds to the amount of point which we took, so, the first attempt would be just do divide that correlation sum to N, let’s see what will happen:

clf
clear
#input data variables
step = 1/1024;
tstop = 1;
t = 0:step:tstop;
y = 2 + 5*cos(2*pi*t*2) + 3*sin(2*pi*t*3) + 2*cos(2*pi*t*4);
#plot(t,y,'b', 'linewidth',4);
grid on;
hold  on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;

tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
#plot(tsample,y_sampled,'r','marker','o','linewidth',1,'markersize',20);

h=legend('original              ','sampled                 ');
set(h,'fontsize',14,'fontname','FreeSans','fontweight','normal');
set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#xlabel('t, s');
#ylabel('Amplitude, V');

for i=0:1:(N-1)
  ytry=e.^(-j*2*pi*i*Fs/N*tsample);
  y_product = (y_sampled.*ytry)/N;
 
  correlation_real(i+1) = sum(real(y_product));
  correlation_imag(i+1) = sum(imag(y_product));
  correlation_abs(i+1) = sqrt(correlation_real(i+1)^2 + correlation_imag(i+1)^2);
  #plot(i,correlation_abs(i+1),'marker','o');
  #hold on;  
  #grid on;
  
  
endfor

xlabel('F, Hz');
ylabel('Amplitude, V');

bar(f_spectrum,correlation_abs,0.1)

OK, a couple of things from my questions were addressed there, but if you draw enough of attention, it could be seen that we have proper values for DC and for F = 4Hz, while 2 and 3 Hz are half from what they should be. Also, we see another two harmonics after 4Hz – this is because we entered in the second Nyquist zone (>Fs/2) and we got the spectrum folding there (and all entries there are just complex conjugates for entries from the other side). So, for instance, what we see at 6 Hz is just 8-2 Hz and so on. This is also the reason why we see /2 degradation inside of the first Nyquist zone, but DC and Fs/2 actually matched in folding points. The conclusion is the next – we need to multiply by 2 points inside of (DC, Fs/2) but not touch edges. Also, we are not interested in the second half of points for the spectrum crafting.

clf
clear
#input data variables
step = 1/1024;
tstop = 1;
t = 0:step:tstop;
y = 2 + 5*cos(2*pi*t*2) + 3*sin(2*pi*t*3) + 2*cos(2*pi*t*4);
#plot(t,y,'b', 'linewidth',4);
grid on;
hold  on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;

tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
#plot(tsample,y_sampled,'r','marker','o','linewidth',1,'markersize',20);

h=legend('original              ','sampled                 ');
set(h,'fontsize',14,'fontname','FreeSans','fontweight','normal');
set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#xlabel('t, s');
#ylabel('Amplitude, V');

for i=0:1:(N-1)
  ytry=e.^(-j*2*pi*i*Fs/N*tsample);
  if ((i==0)||(i==N/2))
    y_product = (y_sampled.*ytry)/N;
  else 
    y_product = (y_sampled.*ytry)/(0.5*N);
  endif  
  correlation_real(i+1) = sum(real(y_product));
  correlation_imag(i+1) = sum(imag(y_product));
  correlation_abs(i+1) = sqrt(correlation_real(i+1)^2 + correlation_imag(i+1)^2);
  #plot(i,correlation_abs(i+1),'marker','o');
  #hold on;  
  #grid on;
  
  
endfor


subplot(3,1,1);
title("mag");
bar(f_spectrum,correlation_abs,0.1)
subplot(3,1,2);
title("real")
bar(f_spectrum,correlation_real,0.1)
subplot(3,1,3)
title("imag")
bar(f_spectrum,correlation_imag,0.1)

With that code magnitude values matched truly across dc – Fs/2, also I have plotted real and imag values to show that the phase information available as well. Also, it is clearly seen that after Fs/2 we have the complex conjugate (take a look on imag plot).

This is great, but there is the one thing still which not quite aligned, if we look on the basic DFT equation:

F[n] = \sum_{k=0}^{N-1}f[k]e^{-j2\pi kn/N} (n=0:N-1)

Octave or matlab allows us to calculate correlation sum just by multiplying waveforms (or matrix) while not all other programming languages can allow that I guess. So, instead, we could “probe” each harmonic (n) by calculation of the correlation sum of the original signal with kn/N

Lets modify the code a little bit, to be compatible with the original DFT equation:

clf;
clear;
grid on;
hold  on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;

tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
#plot(tsample,y_sampled,'r','marker','o','linewidth',1,'markersize',20);

h=legend('original              ','sampled                 ');
set(h,'fontsize',14,'fontname','FreeSans','fontweight','normal');
set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#xlabel('t, s');
#ylabel('Amplitude, V');

for i=0:1:(N-1)
  sum_harm = 0;
  for k=0:1:(N-1)
      sum_harm = sum_harm + y_sampled(k+1)*e^(-j*2*pi*i*k/N);
  endfor
    f(i+1) = sqrt(real(sum_harm)^2+imag(sum_harm)^2);
 
endfor

#normalization procedure

for i=1:1:N
    if ((i==1)||(i==(N/2 + 1)))
      f_norm(i) = f(i)/N;
    else
      f_norm(i) = f(i)*2/N;
    endif
endfor

bar(f_spectrum,f_norm,0.1)
xlabel('f, Hz');
ylabel('Amplitude, V');

And the result is the same:

OK, so DFT is more or less clear I think, there are few notes should be taken into account still:

• Aliasing – that one we have seen already – we have “spectral folding” after Fs/2. To solve that one we can just not use half of the data
• Spectral leakage – in an example I used all harmonics were lying right in the spectral resolution (which is 1Hz for the example), but what if there will be some frequency which has some fractional part

So, lets assume we have 3.3Hz frequency with an amplitude 3 and Fs is equal to 8Hz, also I did N equal to 64, so the resolution is 0.125Hz and our harmonic doesn’t fall in that resolution.

It could be seen that amplitude spreads around somewhat of 3.3Hz and across different fft seeds.

To make this problem less hurting for the resulting spectrum, there is thing called window, which has direct impact on the spectrum by, basically, applying some extra transfer function to the sampled signal being processed.

Let say we have 64 samples, one of the most often used windows is the Hanning one, it has the next function:

w(n) = 0.5(1-cos(2\pi \frac{n}{N}))

It looks like this:

Little bit of a practice:

N = 64;
Fs = 8;
ts = 1/Fs;
#f_spectrum = 0:Fs/N:(Fs*(N-1))/N;
f_spectrum = 0:Fs/N:Fs/2;

tsample = 0:ts:(N-1)*ts;
nsample = 0:1:N;
y_sampled = 2 + 3*cos(2*pi*tsample*3.3);

for i=1:1:(N)
  ywindowed(i) = 2+ 0.5*(1-cos(2*pi*i/N))*(y_sampled(i)-2);
endfor

plot(tsample,y_sampled,'b','linewidth',2,'markersize',10);
hold on;
plot(tsample,ywindowed,'r','marker','o','linewidth',2,'markersize',10);

Here, I just added some post processing to remove dc and then add dc, just to have a nicer plot (which impacts only dc component).

Now, if we apply that Hann window and plot the spectrum we would get:

It is still far from an ideal spectrum, but it is obvious that all extra spectral components which were spreaded across whole spectrum, now concentrated somewhere around 3.3Hz, so the leakage of spectral content away from the correct location is much reduced.

Octave has build in window functions, so instead of code I used, we can use that built-in possibility:

  ywindowed = y_sampled.*hanning(N)';

OK, that one also looks clean more or less I think now. What is left there is the computing power which is used on the calculation of DFT, there are a lot of multiplications made (N*N), more resolution you want more multiplication you have. For instance, if some would want to have a nice spectrum with 65536 points it will result in 4294967296 – quite a significant amount of computations. That’s why FFT was invented which allows reusing already calculated coefficients using decimation in time or decimation in frequency algorithm, but that is out of the scope of this topic and, it is by itself quite a big thing to dig in.

But, I still want to show how to use built-in fft function in Octave, rather than invent own DFT or FFT code:

clf
clear
#input data variables
tstop = 1;
hold  on;
N = 8;
Fs = 8;
ts = 1/Fs;
f_spectrum = 0:Fs/N:(Fs*(N-1))/N;

#define window there
#window  = blackman(N);
window = 1;


tsample = 0:ts:(N-1)*ts;
y_sampled = 2 + 5*cos(2*pi*tsample*2)+3*sin(2*pi*tsample*3) + 2*cos(2*pi*tsample*4);
yFFT = y_sampled.*window';

set(gca, 'linewidth', 3, 'fontsize', 18,'fontname', 'FreeSans') #modify line and fontsize of x and y axes
#plot(tsample,y_sampled,'b','linewidth',2,'markersize',20);
grid on;


spectrum_complex = fft(yFFT, N)/(0.5*N);
spectrum_complex(1) = spectrum_complex(1) /2;
spectrum_complex(N/2+1) = spectrum_complex(N/2+1) /2;

spectrum_mag = abs(spectrum_complex);
bar(f_spectrum(1:(N/2+1)), spectrum_mag(1:(N/2+1)),0.2);

I am not quite sure if I did it right for edge points, but in that example it worked out.