%  gausdemo.m                   T.L.Szabo    Version 1          9/12/02
% this program compares DFTs (FFT) and continuous FTs of Gaussian waveforms 
% for Diagnostic Ultrasound Imaging: Inside Out
% 
close all;
clear all;
delt0=0.1;            %  time sample increment
nt=128;               %  number of time points
t=delt0*(0:nt-1);     %  time array (us)
fc=5.0;               %  center frequency of pulse
a=0.5;                %  exponential Gaussian decay factor
fd0=1./(nt*delt0);    %  frequency increment
f1=fd0*[-nt/2:1.0:nt/2-1];         % sets final frequency range
% 
% sets unshifted  frequency range
fun1=zeros(1,nt);
funa=[0.0:nt/2-1];
funb=[-nt/2:-1];
fun1(1:nt/2)=fd0*funa;         
fun1(nt/2+1:nt)=fd0*funb;         % sets unshifted  frequency range
%
%
% create Gaussian pulse in time
tdelay=nt*delt0/2.0;                   %  shift to middle of range
td=tdelay*ones(1,nt);                  %  create delay array
tt=t-td;
pulse3=exp(-a*tt.*tt);                 %  calculate Gaussian time signal
% 
figure;
plot(t,pulse3);
xlabel('Time (us)');
Title('Gaussian signal');
%
%  calculate DFT spectrum magnitude
spect3=abs(fftshift(fft(delt0*pulse3)));  %take FFT,mult. by delt0,shuffle
figure;
plot(f1,spect3);
Title('DFT Gaussian spectral magnitude');
xlabel('Frequency (MHz)');
%  calculate DFT spectral phase
spect3cmplx=fft(delt0*pulse3);
spect3cor=spect3cmplx.*exp(i*2*pi*tdelay*fun1);
spect3shuf=fftshift(spect3cor);
ang3=unwrap(angle(spect3shuf));
figure;
plot(f1,ang3);
Title('DFT Gaussian spectral phase');
xlabel('Frequency (MHz)');
%
%  calculate continuous Gaussian spectrum obtained analytically
s3=2.506*exp(-19.719*f1.*f1);
figure;
plot(f1,abs(s3));
Title('continuous Gaussian spectrum');
xlabel('Frequency (MHz)');
angle3=unwrap(angle(s3));
figure;
plot(f1,angle3);
Title('continuous Gaussian spectral phase');
xlabel('Frequency (MHz)');
% the end