% Problem 1: Real-Time Fast Convolution using Overlap-Add Technique
%
% The idea is instead of performing convolution in the time domain, which is slow and processor intensive,
% we can take the Fast Fourier Transform of both the input X[n] and the filter H[n] perform convolution
% in the frequency domain, which is simply multiplication, take the inverse fast fourier transform of
% the result and be finished.

% Overlap Add Method (see wikipedia)
% M = 7 - filter size
% N = 10 - sample size
% L = N + M - 1 = 16 - FFT matrix size
% overlaps by (M - 1)  6 samples for each DFT evaluation


% h[n] from fdatool
% 0.0567    0.1560    0.2394    0.2719    0.2394    0.1560    0.0567


h(1) = 0.0567;
h(2) = 0.1560;
h(3) = 0.2394;
h(4) = 0.2719;
h(5) = 0.2394;
h(6) = 0.1560;
h(7) = 0.0567;

% get first 30 samples of function with sample rate of 8khz
pi = 3.141592;
for i=1:30
    t = i * (2 * pi / 8000);
    x(i) = sin(2 * pi * 800 * t) + sin(2 * pi * 1400 * t);
end

% this is what the answer should look like
x_conv = conv(h,x);

subplot(2,1,1), stem(x);
xlabel 'X[n] input'
subplot(2,1,2), stem(h);
xlabel 'H[n] filter'

pause


% zero pad to 36
x(31) = 0;
x(32) = 0;
x(33) = 0;
x(34) = 0;
x(35) = 0;
x(36) = 0;

% zero pad impulse response
for i=8:16
    h(i) = 0;
end



% Get 3 buffers of ten samples
for(i = 1:30)
    if (i <= 10)
        x1(i) = x(i);
    end
    if (i > 10 && i <= 20)
        x2(i-10) = x(i);
    end
    if (i > 20 && i <= 30)
        x3(i-20) = x(i);
    end
end

% zero pad for 16 point DFT
for(i = 11:16)
    x1(i) = 0;
    x2(i) = 0;
    x3(i) = 0;
end




% calculate omega
N=16;
for i = 1:power(N - 1, 2) + 1
    omega(i) = power( exp(-2 * pi * sqrt(-1) / N), i - 1);
end

% populate DFT matrix with omega -- 9x9 matrix for N=10
for(i=1:N)
    for(j=1:N)
        matrix(i,j) = omega( (i-1) * (j-1) + 1 );
    end
end




%y1 = x1 * matrix;
%y2 = x2 * matrix;
%y3 = x3 * matrix;
y1 = fft(x1);
y2 = fft(x2);
y3 = fft(x3);

% So we have the DFT of X[n], do the same for H[n]
%impulse = h * matrix;
impulse = fft(h);


% perform filter in frequency domain
for(i = 1:16)
    y1_filtered(i) = y1(i) * impulse(i);
    y2_filtered(i) = y2(i) * impulse(i);
    y3_filtered(i) = y3(i) * impulse(i);
end


% Inverse is just the complex conjugate of each matrix entry eg: positive sqrt(-1) becomes negative sqrt(-1) and vice versa
inverse_matrix = conj(matrix);

x1_filtered = ifft(y1_filtered);
x2_filtered = ifft(y2_filtered);
x3_filtered = ifft(y3_filtered);

%x1_filtered = y1_filtered * inverse_matrix;
%x2_filtered = y2_filtered * inverse_matrix;
%x3_filtered = y3_filtered * inverse_matrix;

% zero pad
x1 = [x1_filtered zeros(1, 20)]
x2 = [zeros(1, 10) x2_filtered zeros(1, 10)]
x3 = [zeros(1, 20) x3_filtered]

% recombine into single output with overlap adding
for(i = 1:36)
        x_filtered(i) = x1(i) + x2(i) + x3(i);
end


subplot(2,1,1), stem(x_filtered);
xlabel 'filtered time domain result'
subplot(2,1,2), stem(x_conv);
xlabel 'direct convolution'
clear all
pause