# Discrete Fourier Transform
## Artemii Bykov, DS-01

###### tags: `DSP`, `SciLab`, `FFT`, `DFT`, `Spectral leakage`, `Spectral interpolation`

[ToC]

### Task 1
#### DFT Implementation

function [output] = my_dft(data)
    N = length(data);
    for k = 0:N-1
        s = 0;
        for n = 0:N-1
            s = s + data(n+1) * exp(-2 * %i * %pi * n * k / N);
        end
        output(k+1) = s;
    end
endfunction

#### FFT Implementation
Using FFT, we reduce time complexity from $O(N^2)$ to $O(N\log N)$. I have been implemented the most common versions of FFT using a recursive approach.
* Base case of recursion - length of vector equal to 1, just return the vector.
* Otherwise:
* We split the vector into two parts, in one part we have the even indexes, in another the odd indexes
* Compute new additive exponential part for element
* To each current value in vector add multiplication of it is value and the corresponding exponential part
* Concatenate parts and returns
Let is see on the graph of errors
* First error formula: $E_1(x_i) = |Y_i - \hat{Y_i}|$
* Second error formula: $E_2(x_i) = ||Y_i| - |\hat{Y_i}||$
* $MSE = 4.885E-29$
We can see that error is very small but it is not zero. The main reason is the floating-point calculation. In-built FFT uses C-code for calculation, our implementation uses Scilab (in fact Java floating-point calculation)

function [output] = my_fft(data)
    output = conj(fftReq(data)');
endfunction

function [output] = fftReq(data)
    N = length(data);
    if N == 1 then 
        output = data;
    else
        data_even = fftReq(data(1:2:N));
        data_odd = fftReq(data(2:2:N));
        fact = exp(-2 * %i * %pi / N * (0:N-1));
        left = [data_even + fact(1:N/2) .* data_odd];
        right = [data_even + fact(N/2+1:N) .* data_odd];
        output = [left, right];
    end
endfunction

### Task 2
As we have in the tutorial, not integer frequency leads to spectral leakage. We have 3 sinusoids components:
* The same:
* $T$ from $0.001$ to $1$ with frequnecy sampling $1/100$, i.e. length is eqauls to $100$
* The difference:
* The first frequency equals $20.5$
* The second frequency equals $5$
* The third frequency equals $10.5$
The first and the third sinusoids have the spectral leakage, the second one not, it leads to spectral leakage in resulting signal (sum of sinusoids)

**If we linearly increase amplitudes of sinusoids in the time domain, we will have linearly increasing of the magnitude of frequency in the frequency domain and vice versa**

clc;
clear;

fontsize = 3

t = 0.0001 : 1/100:1;

fm1 = 20.5;
x1 = sin(2 * %pi * fm1 * t);

fm2 = 5;
x2 = sin(2 * %pi * fm2 * t);

fm3 = 10.5;
x3 = sin(2 * %pi * fm3 * t);

result_x = x1 + x2 + x3

spectral_leakage = abs(fft(result_x))

result_N = length(result_x)
k = 0:result_N-1;
f = k / result_N

figure(1)
subplot(211)
plot(t, result_x)
xlabel('time', 'fontsize', fontsize), ylabel( 'x(t)', 'fontsize', fontsize)
title('Signal with three sinusoidal components', 'fontsize', fontsize)
subplot(212)
plot(f, spectral_leakage)
title('Spectral leakage', 'fontsize', fontsize)
xlabel('Freq in Hz', 'fontsize', fontsize), ylabel('Mag', 'fontsize', fontsize);

### Task 3
As we have, spectral leakage occurs when a non-integer number of periods of a signal is sent to the DFT. Spectral leakage lets a single-tone signal be spread among several frequencies after the DFT operation. This makes it hard to find the actual frequency of the signal.
Trying to find out the actual frequency of the signal, we can use the spectral interpolation. Zero padding creates more bins for frequencies and more accurate (ideal) interpolation. And we can see spectral leakage was almost significantly decreased.
**Interpolation helps us to get more accurate result of the frequency of the sampled signal:**

clc;
clear;

fontsize = 3

t = 0.0001 : 1/100:1;
t_padded = 0.0001 : 1/100:1/100*128;

fm1 = 20.5;
x1 = sin(2 * %pi * fm1 * t);

fm2 = 5;
x2 = sin(2 * %pi * fm2 * t);

fm3 = 10.5;
x3 = sin(2 * %pi * fm3 * t);

result_x = (x1 + x2 + x3)
padded_result_x = resize_matrix(result_x, 1, 128)


spectral_leakage = abs(fft(result_x))
padded_spectral_leakage = abs(fft(padded_result_x))

result_N = length(result_x)
k = 0:result_N-1;
f = k 

result_N = length(padded_result_x)
k = 0:result_N-1;
padded_f = k / result_N * 100


figure(1)
subplot(221)
plot(t, result_x)
xlabel('time', 'fontsize', fontsize), ylabel( 'x(t)', 'fontsize', fontsize)
title('Signal with three sinusoidal components', 'fontsize', fontsize)
subplot(222)
plot(f, spectral_leakage)
title('Spectral leakage', 'fontsize', fontsize)
xlabel('Freq in Hz', 'fontsize', fontsize), ylabel('Mag', 'fontsize', fontsize);

subplot(223)
plot(t_padded, padded_result_x)
xlabel('time', 'fontsize', fontsize), ylabel( 'x(t)', 'fontsize', fontsize)
title('Signal with three sinusoidal components  (padded)', 'fontsize', fontsize)
subplot(224)
plot(padded_f, padded_spectral_leakage)
title('Spectral leakage (padded)', 'fontsize', fontsize)
xlabel('Freq in Hz', 'fontsize', fontsize), ylabel('Mag', 'fontsize', fontsize);

### Task 4
After compuring FFT, we can see that frequency for the 1st coside is not valid, because sampling rate ($f_s = 200 S/s$) almost eqauls to frequency ($f_1 = 190 Hz$).
And as a consequence, we can observe frequency on the 1st picture equals to $10$, but in fact, it should be equals to $190$.
Now we can use *Nyquist sampling theorem*, which says us, that the sampling
frequency should be at least twice the highest frequency contained in the signal: $f_s \geqslant 2f_c$
I set sampling rate ($f_s = 380 S/s$) for the 1st signal and now we case see frequency equals to 190 (that is true):

clc;
clear;

fontsize = 3


fm1 = 190;
fs1 = 380;
t1 = 0 : 1/fs1 : 0.5;
x1 = 0.5 * cos(2 * %pi * fm1 * t1);
freq1 = t1 ./ t1(length(t1)) * fs1;
fft1 = abs(fft(x1));

fm2 = 10;
fs2 = 200;
t2 = 0 : 1/fs2 : 0.5;
x2 = 2 * cos(2 * %pi * fm2 * t2);
freq2 = t2 ./ t2(length(t2)) * fs2;
fft2 = abs(fft(x2));


figure(1),
subplot(221),
plot2d(t1, x1),
xlabel('Time T', 'fontsize', fontsize),
ylabel('cos(t)', 'fontsize', fontsize),
title('T = 0.5, F = 190Hz, A = 0.5, Fs = 380 S/s', 'fontsize', fontsize),
subplot(223),
plot2d(t2, x2),
xlabel('Time T', 'fontsize', fontsize),
ylabel('cos(t)', 'fontsize', fontsize);
title('T = 0.5, F = 10Hz, A = 2, Fs = 200 S/s', 'fontsize', fontsize);
subplot(222),
plot2d(freq1, fft1),
xlabel('frequency f', 'fontsize', fontsize),
ylabel('magnitude(f)', 'fontsize', fontsize),
title('FFF from T = 0.5, F = 190Hz, A = 0.5, Fs = 380 S/s', 'fontsize', fontsize),
subplot(224),
plot2d(freq2, fft2),
xlabel('frequency f', 'fontsize', fontsize),
ylabel('magnitude(f)', 'fontsize', fontsize),
title('FFT from T = 0.5, F = 10Hz, A = 2, Fs = 200 S/s', 'fontsize', fontsize),

### Bonus task
$$
\frac{1}{N}\sum_{n=0}^{N-1}(e^{-j2\pi k\frac{n}{N}})^*e^{-j2\pi \hat{k}\frac{n}{N}} = \\
= \frac{1}{N}\sum_{n=0}^{N-1}e^{j2\pi \frac{n}{N}(k-\hat{k})} \\
k = \hat{k} \implies e^{j2\pi \frac{n}{N}(k-\hat{k})} = e^0 = 1 \implies \\
\frac{1}{N}\sum_{n=0}^{N-1}1 = \frac{N}{N} = 1 \\
k \neq \hat{k} \implies \frac{1}{N}\sum_{n=0}^{N-1}e^{j2\pi \frac{n}{N}(k-\hat{k})} = \frac{1-e^{j2\pi(k-\hat{k})}}{1-e^{\frac{j2\pi}{N}(k-\hat{k})}} \\
e^{j2\pi} = cos(2\pi) + jsin(2\pi) = 1 = e^0 \\
e^{j2\pi (k-\hat{k})} = e^{0(k-\hat{k})} = 1 \implies \\
\frac{1-e^{j2\pi(k-\hat{k})}}{1-e^{\frac{j2\pi}{N}(k-\hat{k})}} = \frac{1-1}{1-e^{\frac{j2\pi}{N}(k-\hat{k})}} = 0
$$