function cmpcnv(w,sigma);
%This worksheet demonstrates the smoothing effects of convolution and shows
%how to use convolution as a way to approximate derivatives.
%
%From the Command Window use the call cmpcnv(w,sigma)
%
% w is the width of the smoothing kernel. It should be positive and less
% than 1.
% sigma is the amplitude of the noise. .1 is already a lot of noise.
%
%We use 2^k samples
k=10
%sigma is an optional argument giving the noise level
if (nargin < 2) sigma =0; end
%
%Here we choose a function to smooth by uncommenting one:
%
%f = @(x) (abs(x)<1);
%f = @(x) ((cos(pi*x/2)).^2).*(abs(x)<1);
f = @(x) (cos(4*pi*x)>0).*(abs(x)<1);
%
%Next we choose a function to average with:
%
%This does averaging over a rectangle
g= @(t,x) (abs(x)<t)./t;
%This smooths with a twice differentiable function
g=@(t,x) ((t^2-x.^2).^2).*(abs(x)<t);
%This computes the derivative of f*g, where g is the previous function.
g1d= @(t,x) -4*x.*(t^2-x.^2).*(abs(x)<t);
%
%Grid for sampling
%
x=linspace(-2,2,2^k);
ls2=2^(k-1);
%Grid of twice the length for the convolution
x2=linspace(-4,4,2^(k+1));
%Fix the width for the smoothing function here
gw=@(x) g(w,x);
ig=sum(gw(x))*2^(-k-1);
gw=@(x) g(w,x)/ig;
%For derivative uncomment the following
%gw=@(x) g1d(w,x)/ig;
%
%Sample the function and add the noise
fs=f(x)+sigma*randn(size(f(x)));
figure
subplot(1,3,1)
plot(x,fs)
title(['f+\sigma n, \sigma = ',num2str(sigma)])
subplot(1,3,2)
plot(x,gw(x))
title(['g_w, w = ',num2str(w)])
fg=conv(fs,gw(x));
subplot(1,3,3)
plot(x,fg(ls2:(ls2+2^k-1)))
title('f*g_w')
end

function fg=conv(f,g)
%We use the FFT to compute the convolution of f and g which are sampled on 
%a grid 
%
%We zero pad to double the length to correctly convolve them over the given
%interval.
%
%Fs and Gs are the zero-padded samples 
ls=length(f)
Fs=[f zeros(1,ls)];
Gs=[g zeros(1,ls)];
%Compute the FFTs
FFs=fft(Fs);
FGs=fft(Gs);
%Multiply and compute the inverse FFT
FFGs=FFs.*FGs;
fg2=ifft(FFGs);
%Divide by the number of samples to correctly normalize the integrals.
fg=fg2(1:2*ls)/(2*ls);
end