%  AppendixA.m                   T.L.Szabo    Version 2 July 2003          
% this program compares DFTs and continuous FTs of Gaussian waveforms 
% creates appendix A figures in Diagnostic Ultrasound Imaging: Inside Out
close all;
clear all;
delt0=.3;            %  time sample increment
nt=32;                %  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
%
nx=4*nt;                    %  extra points to demonstrate sampling
f4=fd0*[-nx/2:1.0:nx/2-1];  % sets final extra frequency range
% 
% sets desired unshifted  frequency range used with  FFT spectra 
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
%
% sets unshifted extra frequency range for demonstration
fun4=zeros(1,nx);
fun4a=[0.0:nx/2-1];
fun4b=[-nx/2:-1];
fun4(1:nx/2)=fd0*fun4a;         
fun4(nx/2+1:nx)=fd0*fun4b;         % sets unshifted extra 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
pulse3a=exp(-a*tt.*tt);                 %  calculate Gaussian time signal% 
figure;    %    Fig A1a
plot(t,pulse3);
xlabel('Time (us)');
Title('Gaussian signal');
%
figure;    %    Fig A1c
yy=diff(pulse3);
plot(yy);
xlabel('number of samples');
Title('differentiated Gaussian signal');
%
figure;    %   Fig A1b
stem(t,pulse3);
xlabel('Time (us)');
Title('Sampled Gaussian signal');
% 
% create sampled Gaussian pulse in time window with extra points
t4=delt0*(0:nx-1);     %  time array (us)
tdelay4=nx*delt0/2.0;                   %  shift to middle of range
td4=tdelay4*ones(1,nx);                  %  create delay array
tt4=t4-td4;
pulse4=exp(-a*tt4.*tt4);  
%
%  calculate extended sampled Gaussian time signal
figure;    
stem(t4,pulse4);
xlabel('Time (us)');
Title('Sampled extended Gaussian signal');
%
%  create repeated sampled time waveform   Fig A3a
pulse2(1:nt/2)=pulse3(nt/2+1:nt);
pulse2(nt/2+1:3*nt/2)=pulse3;
pulse2(3*nt/2+1:5*nt/2)=pulse3;
pulse2(5*nt/2+1:7*nt/2)=pulse3;
pulse2(7*nt/2+1:4*nt)=pulse3(1:nt/2);
figure;              %  FigA3a
stem(t4(nt/2+1:7*nt/2),pulse2(nt/2+1:7*nt/2));
%plot(f4,spect2);
Title('Repeated sampled time waveform');
xlabel('Time (usec)');
%
%  calculate DFT spectrum magnitude Fig A2a
spect3=abs(fftshift(fft(delt0*pulse3)));  %take FFT,mult. by delt0,shuffle
figure;         
plot(f1,spect3);
Title('DFT Gaussian spectral magnitude');
xlabel('Frequency (MHz)');
%
figure;         % Fig A6b
stem(f1,spect3);
Title('DFT Gaussian spectral magnitude');
xlabel('Frequency (MHz)');%  calculate spectral phase
spect3cmplx=fft(delt0*pulse3);
spect3cor=spect3cmplx.*exp(i*2*pi*tdelay*fun1); % correct for delay
spect3shuf=fftshift(spect3cor);
ang3=unwrap(angle(spect3shuf));
figure;  % phase for Fig A6b
plot(f1,ang3);
Title('DFT Gaussian spectral phase');
xlabel('Frequency (MHz)');
%
%  calculate unshifted DFT spectrum magnitude
spect7=abs(fft(delt0*pulse3a));  %take FFT,mult. by delt0
%
figure;
spect6=abs(fftshift(fft(delt0*pulse3)));  %take FFT,mult. by delt0,shuffle
stem(spect7);%
Title('Unshifted DFT Gaussian spectral magnitude');
xlabel('Frequency (MHz)');
%  calculate extended DFT spectrum magnitude
spect4=abs(fftshift(fft(delt0*pulse4)));  %take FFT,mult. by delt0,shuffle
figure;
plot(f4,spect4);
Title('Extended DFT Gaussian spectral magnitude');
xlabel('Frequency (MHz)');
%  calculate spectral phase
%
spect4cmplx=fft(delt0*pulse4);
spect4cor=spect4cmplx.*exp(i*2*pi*tdelay4*fun4);  %correct for delay
spect4shuf=fftshift(spect4cor);
ang4=unwrap(angle(spect4shuf));
figure;
plot(f4,ang4);
Title('Extended DFT Gaussian spectral phase');
xlabel('Frequency (MHz)');
%
% simulate continuous spectrum of a sampled waveform
spect2(1:nt/2)=spect3(nt/2+1:nt);
spect2(nt/2+1:3*nt/2)=spect3;
spect2(3*nt/2+1:5*nt/2)=spect3;
spect2(5*nt/2+1:7*nt/2)=spect3;
spect2(7*nt/2+1:4*nt)=spect3(1:nt/2);
figure;              %  FigA2B
plot(f4(nt/2+1:7*nt/2),spect2(nt/2+1:7*nt/2));
Title('Spectral magnitude of sampled Gaussian');
xlabel('Frequency (MHz)');
%
figure;              %  FigA3b
stem(f4(nt/2+1:7*nt/2),spect2(nt/2+1:7*nt/2));
Title('Spectral magnitude of repeated sampled Gaussian');
xlabel('Frequency (MHz)');
%
%   calculate time delay - corrected pulse
ttd=delt0*(-nt/2:1.0:nt/2-1);     %  time array (us)
pulsed=exp(-a*ttd.*ttd);                 %  calculate Gaussian time signal
%pulseshuf=ifftshift(pulsed);%
figure;    %    Fig A5a
plot(ttd,pulsed);
%plot(ttd,pulseshuf);
xlabel('Time (us)');
Title('Undelayed Gaussian signal');
%   calculate analytic signal
spect3cmplx(nt/2+1:nt)=zeros(1,nt/2);
pulsea=ifft(spect3cmplx)*2/delt0;
pulseh=hilbert(pulsed);
ttd=delt0*(-nt/2:1.0:nt/2-1);     %  time array (us)
pulseshuf=ifftshift(pulsed);%
figure;    %    Fig A7
plot(ttd,real(pulseh),'k-',ttd,imag(pulseh),'k:',ttd,abs(pulseh),'k-.');
%plot(ttd,pulsea);
xlabel('Time (us)');
Title('Undelayed analytic Gaussian signal');
legend('real','imaginary','envelope');
%
%  calculate continuous Gaussian spectrum
s3=2.506*exp(-19.719*f1.*f1);
figure;  % compare continuous spectrum with DFT version fig A2A
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
