% eda09_09
% amplitude spectrum of Normal smoothing filter

% standard set up of time and frequency
N=256;
Dt=0.01;
tmax=Dt*(N-1); 
t=Dt*[0:N-1]'; 
fmax=1/(2.0*Dt);
df=fmax/(N/2);
f=df*[0:N/2,-N/2+1:-1]'; 
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]';
Nw=Nf;

% evaluate fourier transform of filter, g, of length L
L1=3;

% compute variance of filter (as if it were distribution)
g1 = ones(L1,1);
g1 = g1/sum(g1);
mg1 = sum( [1:L1]' .* g1 ) ;
vg1 = sum( (([1:L1]'-mg1).^2) .* g1);
sg1 = Dt*sqrt(vg1);

% fourier transform of gausian, using analytic formula
sigma1=1/sg1;
s1 = exp(-w(1:Nf).*w(1:Nf)/(2*sigma1^2));
sa1 = abs(s1);
disp(sprintf('length %d samples (%f seconds) sigma-t samples %f (%f seconds)', L1, Dt*L1, sg1/Dt, sg1 ));

figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1] );
plot( f(1:Nf), sa1, 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('|s| for L=3');

% evaluate fourier transform of filter of length L
L2=21;

% compute variance of filter (as if it were distribution)
g2 = ones(L2,1);
g2 = g2/sum(g2);
mg2 = sum( [1:L2]' .* g2 ) ;
vg2 = sum( (([1:L2]'-mg2).^2) .* g2);
sg2 = Dt*sqrt(vg2);
disp(sprintf('length %d samples (%f seconds) sigma-t samples %f (%f seconds)', L2, Dt*L2, sg2/Dt, sg2 ));

% fourier transform of gausian, using analytic formula
sigma2=1/sg2;
s2 = exp(-w(1:Nf).*w(1:Nf)/(2*sigma2^2));
sa2 = abs(s2);

subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, fmax, 0, 1] );
plot( f(1:Nf), sa2, 'k-', 'LineWidth', 2);
xlabel('frequency, Hz');
ylabel('|s| for L=21');