%  blockdemo.m       T.L.Szabo            July 2004
% calculates the output of  a filter for an input waveform & its spectrum
% this program creates figures for waveforms and their spectra for
% Chapter 2 in Diagnostic Ultrasound Imaging: Inside Out% 
close all;
clear all;
delt0=0.01;           %  time sample increment
nt=512;               % number of time points
t=delt0*(1:nt);       %  time array (us)
fc=5.0;               %  center frequency of pulse
a=20.0;               %exponential Gaussian decay factor
d=50;               %exponential Gaussian decay factor
%
%B
% 5 MHz Gaussian pulses: input (pulse1) & filter (pulse2)
tdelay=1.27;
td=tdelay*ones(1,nt);
tt=t-td;
tt2=t-2*td;
pulse1=exp(-d*tt.*tt).*cos(2*pi*fc*tt);
pulse2=exp(-a*tt2.*tt2).*cos(2*pi*fc*tt2);
%
% calculate time waveforms (input & filter
figure;
plot(t,pulse1,'r',t,pulse2,'k');
xlabel('Time (us)');
ylabel('Amplitude');
legend('input','filter');
Title('input and filter pulses');
%
%  calculate shifted and unshifted frequency arrays
fd0=1./(nt*delt0);     %  frequency increment
f1=fd0*[-nt/2:1.0:nt/2-1];            % sets frequency range
% sets unshifted  frequency range corresponding to DFT calculation
fun1=zeros(1,nt);                % sets up initial dummy array of zeros
funa=[0.0:nt/2-1];               % positive frequency indices
funb=[-nt/2:-1];                 % negative frequency indices
fun1(1:nt/2)=fd0*funa;           % equates positive frequency segment
fun1(nt/2+1:nt)=fd0*funb;        % equates positive frequency segmentspect1un=zeros(1,nt);                 % initialize spectrum
%
%  calculate spectra
spect1un=zeros(1,nt);                 % initialize spectrum
spect1cmplx=delt0*fft(pulse1);        % take fft of pulse1spect1un=zeros(1,nt);                 % initialize spectrum
spect1un(1:nt/2)=spect1cmplx(1:nt/2); %set negative freq. spectra to =0spect3
spect1=fftshift(spect1un);            %  shuffle data to standard freq.'s
spect2un=zeros(1,nt);                 % initialize spectrum
spect2cmplx=delt0*fft(pulse2);        % take fft of pulse2
spect2un(1:nt/2)=spect2cmplx(1:nt/2); %set negative freq. spectra to =0
spect2=fftshift(spect2un);            %  shuffle data to standard freq.'s
spect3cmplx=spect1cmplx.*spect2cmplx; % multiply spectra to find output
%spect3cor=spect3cmplx.*exp(i*2*pi*tdelay*fun1); % remove filter tdelay shift
spect3cor=spect3cmplx;                %  include true filter delay shift
spect3un=zeros(1,nt);                 % initialize spectrum
spect3un(1:nt/2)=spect3cor(1:nt/2);   %set negative freq. spectra to =0spect3
spect3=fftshift(spect3un);            %  shuffle data to standard freq.'s
% plot spectra
figure;
plot(f1,abs(spect1),'r--',f1,abs(spect2),'k',f1,abs(spect3),'m');
Title('pulse spectral magnitudes'); xlabel('Frequency (MHz)');
ylabel('Amplitude');legend('input','filter','output');
axis([-10.  10.0 0.0 0.5]);
%
%  calculate inverse fft to get analytic signal
pulse1a=2*ifft(spect1un)/delt0;
pulse1env=abs(pulse1a);
pulse2a=2*ifft(spect2un)/delt0;
pulse2env=abs(pulse2a);
pulse3a=2*ifft(spect3un)/delt0;
pulse3env=abs(pulse3a);
%
figure;
plot(t,real(pulse1a),'r--',t,real(pulse2a),'k',t,real(pulse3a),'m');
hold;
plot(t,pulse1env,'r--',t,pulse2env,'k',t,pulse3env,'m');
xlabel('Time (us)');ylabel('Amplitude');
legend('input','filter','output');Title('real pulses & their envelopes');
hold;
% the end
