% eda06_11
% integral using fft example

clear all;

% standard set up
N=256;
Dt=1.0;
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;

% exemplary function
t1 = tmax/4;
t2 = 3*tmax/4;
s = 25;
d = exp( -((t-t1).^2) / (2*s^2) ) - exp( -((t-t2).^2) / (2*s^2) );

% derivative using Riemann sum
integral1 = Dt*cumsum(d);

% derivative using fft
integral=ifft(-i*fft(d).*[0,1./w(2:N)']');

figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -1, 1] );
plot(t,d,'k-','LineWidth',2);
xlabel('time, t');
ylabel('d(t)');
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -100, 150] );
plot(t,integral1,'k-','LineWidth',2);
xlabel('time, t');
ylabel('integral');
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, tmax, -100, 150] );
plot(t,real(integral),'k-','LineWidth',2);
xlabel('time, t');
ylabel('integral');








