% program to perform beamformation for a phased array
% beamforming scheme: RF interpolation
% output format: baseband demodulated IQ data

% define data file parameters
infile='acuson17c.dat';
Nsample=2048;
Nchan=128;
Nskip=0;
Nrecord=1;

% transducer parameters
pitch=0.22;
soundv=1.48;
timeoffset=29.448;
f0=3.5;
fad=13.8889;

% compute delay data
display('compute delay data');
Nrange=256;
tfnumber=2;
rfnumber=1.5;
tfocus=60;
tapersize=tfocus/tfnumber/2;
apersize=Nchan*pitch;
lambda=soundv/f0;
dsin=lambda/apersize/2;
Nbeam=round(sqrt(2)/dsin);
drange=soundv/fad/2;
rangeoffset=timeoffset*soundv/2;
upsamp=8;
beamsum=zeros(Nbeam,Nsample);
%%%%%%%%%%%%%%
for iloop=1:8
%%%%%%%%%%%%%%
iloop 
fid=fopen(infile,'r');
[a,count]=fread(fid,Nskip*Nsample,'integer*2');
Nstart=(iloop-1)*Nrange+1;
Nend=Nstart+Nrange-1;
tdelay=-9999*ones(Nchan,Nbeam);
rdelay=-9999*ones(Nchan,Nbeam,Nrange);
for i=1:Nchan
   x=(i-(Nchan+1)/2)*pitch;
   ax=abs(x);
   tdelayf=x*x/2/tfocus/soundv;
	for j=1:Nbeam
   	sint=dsin*(j-(Nbeam+1)/2);
   	cost2=1-sint*sint;
      cost=sqrt(cost2);
      if ax<=tapersize,
         temp=tdelayf*cost2-x*sint/soundv;
         tdelay(i,j)=round(temp*fad*upsamp);
      end
   	rangeon=rfnumber*2*abs(x)*cost;
   	for k=1:Nrange
         range=rangeoffset+(k+Nstart-2)*drange;
         if range>=rangeon,
            temp=(x*x*cost2/2/range-x*sint)/soundv;
            rdelay(i,j,k)=round(temp*fad*upsamp);
         end
      end
   end
end

% beamformation master loop
display('beam forming');
beamout=zeros(Nbeam,Nrange);
delaysum=beamout;
alpha1=beamout;
alpha2=beamout;
i1=beamout;
for i=1:Nchan
   i
   for j=1:Nchan
      [a,count]=fread(fid,Nrecord*Nsample,'integer*2');
      a=[a;0;0];
      delaysum=tdelay(i,:).'*ones(1,Nrange)+shiftdim(rdelay(j,:,:));
      i1=ones(Nbeam,1)*[1:Nrange]+Nstart-2+floor(delaysum/upsamp);
      alpha2=mod(delaysum,upsamp)/upsamp;
      alpha1=1-alpha2;
      inrange=i1>0 & (i1+1)<Nsample;
		i1=i1.*inrange+(1-inrange)*(Nsample+1);
      beamout=beamout+alpha1.*a(i1)+alpha2.*a(i1+1);
   end
end
beamsum(:,Nstart:Nend)=beamout;
ft=fclose(fid);
%%%%%
end
%%%%%
save out2 beamsum