% Circular Convolution

% NOTE: This demo originally computed the number of floating point operations
%   for each filtering process using "flops".  Matlab 6 no longer supports 
%   this function, so this aspect of the demo was removed.
%   "cputime" doesn't have adequate resolution to substitute for this.

%---------------------------------------------------------------
% Initialize
clear all;
close all;
%---------------------------------------------------------------
% Parameters
M = 256		% window length
N = 1024	% signal length
a = 1./(N-1)
b = 0.0
c = 0.5
x = a.*(0:1:N-1) + b + c.*randn(1,N);
%h = hamming(M)'/M;
whamming=.54-.46*cos(2*pi/(M-1)*(0:M-1));
h=whamming/M;

%---------------------------------------------------------------
% Filtering by direct convolution
%---------------------------------------------------------------
flops(0)
y = conv(x,h);
%C = flops

figure(1)
subplot(3,1,1); plot(h); axis([1 N+M-1 0 2/M])
%ttext = sprintf('Filtering by direct convolution  (C = %d)',C)
ttext = sprintf('Filtering by direct convolution')
title(ttext)
ylabel('h')
subplot(3,1,2); plot(x); axis([1 N+M-1 -2 2])
ylabel('x')
subplot(3,1,3); plot(y); axis([1 N+M-1 0 1])
ylabel('y')
pause

%---------------------------------------------------------------
% Filtering by N pt. FFT
%---------------------------------------------------------------
h(M+1:N) = zeros(1,N-M);

%flops(0)
X = fft(x);
H = fft(h);
Y = H.*X;
y = ifft(Y);
%C = flops

figure(2)
subplot(3,1,1); plot(h); axis([1 N+M-1 0 2/M])
%ttext = sprintf('Filtering by N pt. FFT  (N = %d, C = %d)',N,C)
ttext = sprintf('Filtering by N pt. FFT  (N = %d)',N)
title(ttext)
ylabel('h')
subplot(3,1,2); plot(x); axis([1 N+M-1 -2 2])
ylabel('x')
subplot(3,1,3); plot(real(y)); axis([1 N+M-1 0 1])
ylabel('y')
pause

%---------------------------------------------------------------
% Filtering by N+M-1 pt. FFT
%---------------------------------------------------------------
x(N+1:N+M-1) = zeros(1,M-1);
h(N+1:N+M-1) = zeros(1,M-1);

%flops(0)
%flops
X = fft(x);
H = fft(h);
Y = H.*X;
y = ifft(Y);
%C = flops


figure(3)
subplot(3,1,1); plot(h); axis([1 N+M-1 0 2/M])
%ttext = sprintf('Filtering by N+M-1 pt. FFT  (N+M-1 = %d, C = %d)',N+M-1,C)
ttext = sprintf('Filtering by N+M-1 pt. FFT  (N+M-1 = %d)',N+M-1)
title(ttext)
ylabel('h')
subplot(3,1,2); plot(x); axis([1 N+M-1 -2 2])
ylabel('x')
subplot(3,1,3); plot(real(y)); axis([1 N+M-1 0 1])
ylabel('y')
pause

%---------------------------------------------------------------
% Filtering by N+M pt. FFT
%---------------------------------------------------------------
x(N+1:N+M) = zeros(1,M);
h(N+1:N+M) = zeros(1,M);

%flops(0)
%flops
X = fft(x);
H = fft(h);
Y = H.*X;
y = ifft(Y);
%C = flops


figure(3)
subplot(3,1,1); plot(h); axis([1 N+M-1 0 2/M])
%ttext = sprintf('Filtering by N+M pt. FFT  (N+M = %d, C = %d)',N+M,C)
ttext = sprintf('Filtering by N+M pt. FFT  (N+M = %d)',N+M)
title(ttext)
ylabel('h')
subplot(3,1,2); plot(x); axis([1 N+M-1 -2 2])
ylabel('x')
subplot(3,1,3); plot(real(y)); axis([1 N+M-1 0 1])
ylabel('y')
pause

%---------------------------------------------------------------
% Filtering by L = 2^ceil(log2(N+M-1)) pt. FFT
%---------------------------------------------------------------
L = 2^ceil(log(N+M-1)/log(2))
x(N+M:L) = zeros(1,L-(N+M-1));
h(N+M:L) = zeros(1,L-(N+M-1));

%flops(0)
%flops
X = fft(x);
H = fft(h);
Y = H.*X;
y = ifft(Y);
%C = flops

figure(4)
subplot(3,1,1); plot(h); axis([1 N+M-1 0 2/M])
%ttext = sprintf('Filtering by L = 2^ceil(log2(N+M-1)) pt. FFT  (L = %d, C = %d)',L,C)
ttext = sprintf('Filtering by L = 2^{ceil(log2(N+M-1))} pt. FFT  (L = %d)',L)
title(ttext)
ylabel('h')
subplot(3,1,2); plot(x); axis([1 N+M-1 -2 2])
ylabel('x')
subplot(3,1,3); plot(real(y)); axis([1 N+M-1 0 1])
ylabel('y')