%Pulse echo 1D imaging - simulation
% simulate received signal from an 1 D object defined by variation i acoustic impedance

%19.01.02  Hans Torp
%modified for MEDT8002 13.09.04 HT

%physical constants
c=1540;  %speed of sound (her brukes gjennomsnittsverdi for mykt humant vev)

fs=160e6;
dt=1/fs;
dz=dt*c/2;% depth increment
zmax=0.04;%max depth
z=0:dz:zmax;%depth axis
z=z';%make z into a vector

t=2/c*z;%time axis
tmax=t(end);

% calculate the total pulse echo response h = pel x pxd x hd x hxd ; where "x" means convolution
%define transmitted pulse
f0=5.0e6;
lambda=c/f0;
Tp=0.6e-6;%pulse length
tp=0:1/fs:Tp;tp=tp';
pel=0.3*sin(2*pi*f0*tp);% Electrical transmit pulse

%define transduce impulse response
fc=5.0e6;B=3e6;%center frequency and bandwidth of transducer
[bxd,axd]=butter(2,2*[fc-B/2,fc+B/2]/fs);%butterworth bandpass filter
txd=0:dt:2e-6;txd=txd';
impuls=zeros(size(txd));impuls(1)=1; %
hxd=1*filter(bxd,axd,impuls);
%fla12L_imp;%reading impulse respons from file
txd=(1:length(hxd))*dt;
hxd=hxd/max(hxd);
figure(2);
plot(txd,hxd);

pAc=conv(pel,hxd);%transmitted acoustical pulse
hd = [-1;1];% differensiation operator
h=conv(pAc,hd);%differensiation of Ac. impedance
h=conv(h,hxd);% transducer receiver response = transmit response hxd

%object defined by acoustic impedance Z as a function of depth
% Object 1: water/fat interface in depth z1
Zw=1.48; % [kg/m^2/s] acoustic impedance of water
Zb=1.61; % [kg/m^2/s] acoustic impedance of blood
Zm1=1.63; % [kg/m^2/s] acoustic impedance of muscle
Zm2=1.74; % [kg/m^2/s] acoustic impedance of muscle
Zf=1.37; % [kg/m^2/s] acoustic impedance of human fat tissue

%object defined by acoustic impedance Z as a function of depth
obj=2;
Z=ones(size(z))*Zw;
switch obj
    case 1 % water / fat interface
        z1=0.01;
        Z(find(z>z1))=Zf;
    case 2 % water / fat/ water interface
        z1=0.02; dz1=1e-3;%dz1=lambda/4;
        Z(find(z>z1))=Zf;
        Z(find(z>(z1+dz1)))=Zw;
    case 3 % muscle 2 /  muscle 1 / blood / muscle1/muscle2
        z1=0.02;  dz1=0.5e-3;
        z2=0.027; dz2=0.5e-3;
        Z=ones(size(z))*Zm2;
        Z(find(z>z1))=Zm1;
        Z(find(z>z1+dz1))=Zb;
        Z(find(z>z2))=Zm1;
        Z(find(z>z2+dz2))=Zm2;
    case 4 %muscle/fat random mixture between z1 and z2 in fat
        z1=0.01;z2=0.03;
        Z=ones(size(z))*Zf;
        I12=find(z>z1&z<z2);ism2=rand(size(I12));
        Z(I12)=Zm1*(1-ism2)+Zm2*ism2;
       
end;%switch obj

figure(1);
subplot(1,3,1);plot(Z,z);ylabel('depth z [m]');axis ij;axis([1.3,1.8,z(1),z(end)]);

s=conv(h,Z);%received signal
s=s(1:end-length(h)+1);
an=1e-3;
sn=an*randn(size(s));%Gaussian white noise
s=s+sn;%add thermal noise

smax=max(abs(s));
smax=60;
subplot(1,3,2);plot(s,t);ylabel('time [s]');axis([-smax,smax,t(1),t(end)]);axis ij;
amp=abs(hilbert(s));%amplitude of echo signal s
logamp=20*log10(amp);
subplot(1,3,3);
%plot(t,logamp);xlabel('time [s]');axis([t(1),t(end),-30,30]);
image([0,1],z,[100*(1.8-Z),logamp+0]);colormap(gray(40));
zoom yon
