% Fast convolution demo (overlap and save)
% Nov 30, 2015
% for ECE4703
% -------------------------------------------------------
% User parameters below
% -------------------------------------------------------
N = 512;            % FFT size (must be power of 2)
fs = 44100;         % sampling frequency
NN = 5E5;           % number of input samples to generate (should be big)
% -------------------------------------------------------

% generte input
x = randn(NN,1);    % random input

% design filter
Fpass = 3000;            % Passband Frequency
Fstop = 3500;            % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.01;            % Stopband Attenuation
dens  = 20;              % Density Factor
[h_ord, Fo, Ao, W] = firpmord([Fpass, Fstop]/(fs/2), [1 0], [Dpass, Dstop]);
h = firpm(h_ord, Fo, Ao, W, {dens});
h = h(:);
M = length(h);           % number of filter coefficients (h_ord+1)
K = N-(M-1);

% direct convolution
disp('direct convolution with Matlab...');
tic
y1 = conv(x,h); % <- Matlab cheats!
toc

% direct convolution (this is in no way optimal)
xx = [zeros(M-1,1) ; x ; zeros(M-1,1)]; % zero pad
y2 = zeros(NN+M-1,1);
disp('direct convolution with for loop...');
for n=1:NN+M-1,
    for k=1:M,
        y2(n) = y2(n)+h(k)*xx(n+M-k);
    end
end
toc

% overlap and save fast convolution
disp('fast convolution...');
tic
H = fft([h;zeros(N-M,1)]);      % FFT of zero-padded filter coefficients
y3 = zeros(NN+M-1,1);           % allocate space for output
index = 1;                      % starting index for input
z = zeros(M-1,1);               % initial values for input block prepend
while (index+K-1)<=NN,
    xblock = x(index:index+K-1);    % grab input block of length K
    u = [z;xblock];                 % prepend
    U = fft(u);                     % take FFT
    V = U.*H;                       % multiplication in freq domain
    v = ifft(V);                    % take IFFT
    y3(index:index+K-1) = v(M:end); % throw away first M-1 samples to get output
    z = xblock(end-M+2:end);        % compute new input block prepend
    index = index+K;                % update index
end
toc