%The Fourier transform of a pulse is a sinc, which has infinite bandwidth.
%What is the effect on the original time pulse if we band-limit the sinc
%(which must be done in practice in
%order to obey spectrum allocation rules, for example FCC reg's)?
%In this script, we directly compute the inverse Fourier transform of a
%band-limited sinc using Eq.'s 3.38 and 3.39 in Kamen & Heck.
%This script uses the continuous-time and continuous-frequency definitions
%of the forward and inverse transforms (not the DTFT or DFT). It performs
%the integrations numerically using time and frequency resolutions
%specified by the user.

close all; %close all open figure windows
clear; %delete all variables in memory
tau=1; %full width of original time pulse (units are arbitrary,
% but can assume seconds, ms, us, etc. if desired)
t1=1; %compute inverse Fourier transform (IFT) over -t1 to +t1
wb=32*pi; %truncate sinc to -wb < w < +wb
nts=256; %number of timesteps at which to compute IFT
tvals=linspace(-t1,t1,nts); %time values at which to compute IFT
tinc=tvals(2)-tvals(1); %time increment
nfs=4096; %number of frequency steps
w=linspace(-wb,wb,nfs); %frequency samples
winc=w(2)-w(1); %frequency increment
sincfun=(1./w).*sin(w*tau/2); %sinc function = transform of pulse
subplot(1,2,1)
plot(w,sincfun); %show truncated sinc
xlabel('w, rad/s');title('Bandlimited Fourier Transform of Ideal Pulse')
grid on
n=0; %initialize the loop over time
for t=-t1:tinc:t1 %Compute the IFT integral at each time step
n=n+1;
cexp=exp(1i.*w*t); %complex exponential factor
funw=sincfun.*cexp; %multiply factors to get function of w to be integrated
xp(n)=(1/pi)*trapz(w,funw); %perform numerical integration using trapezoidal rule to recover
% approximation of the original ideal pulse

end
subplot(1,2,2)
plot(tvals,abs(xp)); %plot recovered pulse
xlabel('time');title('Recovered Pulse');grid on