% E03_3.m        17oct14   25mar22
% Linear Convolution Via the DFT    
% Overlap-add method 
clear all; clf;

x = [10 2 -5 4 -5 -3 7 8 -3 9 5];       %  first sequence
h = [2 -1 3];             %  second sequence
L = length(x)+length(h)-1;   % Determine the length of the result of convolution

XE = fft(x,L);               % Apply zero padding and compute the DFT
HE = fft(h,L);
%  XE = fft(x);           % WRONG  ...but can be done!
%  HE = fft(h,length(x)); % (outputs shifted conv., wrong transient, 'same' length  
y1 = ifft(XE.*HE);   %  product and IDFT
L1 = length(y1);

y2 = conv(x,h,'full');   %  compare with convolution, try 'same' and 'full' (default)
L2 = length(y2);

figure(1);
subplot(311);   stem(0:L1-1,y1)
xlabel('Time index n');ylabel('Amplitude');title('DFT-based linear convolution')
subplot(312);   stem(0:L2-1,y2)
xlabel('Time index n');ylabel('Amplitude');title('Time-domain linear convolution')
if L1 == L2,
    subplot(313);  stem(0:L1-1,y1-y2)
    xlabel('Time index n');ylabel('Amplitude'); title('Error sequence')
    end
return;
 
%% Overlap-add method for linear filtering
figure(1);          %  N should be large
N = 1e7;  f1 = 1e6;  f2 = 2e6;   % ang.freq.:  2 pi f / N = .2 pi  and  .4 pi  rad/sample
n = 0:N-1;                  %  two sine functions plus noise:
s = cos(2*pi*n*f1/N) + sin(2*pi*n*f2/N);
d = randn(1,N)*3;  % e.g. up to *10
x = s + d;
 h = [1 3 5 7 9 9 7 5 3 1]/50;        % for a short filter "filter" is faster
 h = sinc(-100:.3:100); h=h./sum(h);  % for a long filter "fftfilt" is faster
tic;  y1 = fftfilt(h,x);  toc      % Overlap-add method of FIR filtering using FFT
tic,  y2 = filter(h,1,x);  toc      % (Direct Form II Transposed)
subplot(221);  plot(s(2001:2100));   %  signals f1 and f2
subplot(222);  plot(x(2001:2100));    %  signals and noise
subplot(223);  plot(y1(2001:2100));  %  filtered: ~ only f1
subplot(224);  plot(y1(2001:2100)-y2(2001:2100));   % difference btwn overlap-add and direct filter

figure(2); % take a look at the DFT of the signals
S = abs(fft(s));     X = abs(fft(x));      Y1 = abs(fft(y1));
subplot(311); plot([0:N-1]/N*2*pi,S);  axis([0,2*pi,0,16e4]);  % spectra
subplot(312); plot([0:N-1]/N*2*pi,X);  axis([0,2*pi,0,16e4]);   % 
subplot(313); plot([0:N-1]/N*2*pi,Y1);  axis([0,2*pi,0,16e4]);   % 
%  NB: note DFT takes large values; 1/N scaling is done in the IDFT
