function PwDoppler;
% Pulsed wave signal from moving reflector
c=1540;
v= 0.2;%m/sec
%v= 0.9625;
%v= 1.92;
prf=5e3;%pulse repetition frequency

z0=0.01;%start position
z1=0.013;%end position

complexDemod=0;%=1 to perform complex demodulation

prt=1/prf;
dz=abs(v*prt);
z=z0:dz:z1;
if v>0,z=fliplr(z);end;
zs=0.0125;
figure(1);clf;
pulseEchoPlot(z(1),zs,complexDemod);
pause;
n=1;
figure(2);clf;

for n=1:length(z)
    figure(1);
    iq(n)=pulseEchoPlot(z(n),zs,complexDemod);
    t(n)=n*prt;
    if 1,
        figure(2);
        if complexDemod,
            subplot(2,1,1);plot(iq);hold on;plot(iq(end),'*');hold off;axis('image');axis([-1,1,-1,1]*1.5);grid;
            xlabel('Real part');ylabel('Imaginary part');
            subplot(2,1,2);plot(t,real(iq),t,imag(iq));
        else subplot(2,1,2);plot(t,real(iq));end;
        axis([0,0.015,-1.5,1.5]);pause(0.01);
    end;
end;
figure(2);clf;
if complexDemod,
subplot(2,1,1);plot(iq);axis('image');axis([-1,1,-1,1]*1.5);grid;
xlabel('Real part');ylabel('Imaginary part');
subplot(2,1,2);plot(t,real(iq),t,imag(iq));
else subplot(2,1,2);plot(t,real(iq));end;
axis([0,0.015,-1.5,1.5]);
Ns=2*prf;n=randn(1,Ns)+i*randn(1,Ns);
siq=filter2(fliplr(iq),n);siq=siq/max(siq);
t1=(1:length(siq))/prf;
pause;
if complexDemod,
    subplot(2,1,1);plot(t1,real(siq),t1,imag(siq));
else
    subplot(2,1,1);plot(t1,real(siq));
end;
zoom xon;
sound(real(siq),prf); 


function iq1=pulseEchoPlot(z1,zs,complexDemod);
%Pulse echo 1D imaging - simulation
% z1 is distance to reflector [m]. default=0.01 m
%23.01.01  Hans Torp
% Modified for complex demodulation 17.10.01

if nargin<1,z1=10e-3;end;

c=1540;%speed of sound
fs=100e6;%sampling freq.
tmax=20e-6;
t=0:1/fs:tmax;%time axis
z=c/2*t;dz=z(2)-z(1);%depth axis

%define  transmit pulse
f0=2e6; B=0.5e6;
Tp=3e-6;%pulse length
tp=0:1/fs:Tp;
p=sin(2*pi*f0*tp);%rectangular pulse
p=XmitPulse(fs,f0,B);
obj=0*z;obj(1)=2;%object function;
obj(1+round(z1/dz))=1;%reflector in position z1
s=filter(p,1,obj);%received signal

figno=1;
figure(figno);set(figno,'DoubleBuffer','on');%prevent flickering

sp=hilbert(s);%pre envelope
sh=imag(hilbert(s));%hilbert transformed signal
smix=exp(-i*2*pi*f0*t);%down-mixing signal
iq=sp.*smix;%in phase and quadrature components
[dum,is]=min(abs(z-zs));
iq1=iq(is);
ts=2/c*zs;

myaxis=[t(1),t(end),-2.2,2.2];
subplot(4,1,1);plot(z,obj);xlabel('depth z [m]');axis('tight');ylabel('Object');
subplot(4,1,2);plot(t,s);xlabel('time [s]');axis(myaxis);ylabel('RF-signal');
hold on; plot([ts,ts],[-1,1]);hold off;
if complexDemod,
    subplot(4,1,3);plot(t,sh);xlabel('time [s]');axis(myaxis);ylabel('Hilbert transf.');
    hold on; plot([ts,ts],[-1,1]);hold off;
    subplot(4,1,4);plot(t,real(iq),t,imag(iq));xlabel('time [s]');axis(myaxis);ylabel('IQ signal');
    hold on; plot([ts,ts],[-1,1]);hold off;
end;

pause(0.01);

function pxd=XmitPulse(fs,fc,B,Tp);
%define transmitted pulse
%fs: sampling frequency
%fc: pulse center frequency
%B  pulse bandwidth
% Tp  pulse length, default Tp=3/B;
if nargin<4, Tp=3/B;end;
[bxd,axd]=butter(2,2*[fc-B/2,fc+B/2]/fs);%butterworth bandpass filter
dt=1/fs;
txd=0:dt:3e-6;txd=txd';
impuls=zeros(size(txd));impuls(1)=1; %
pxd=filter(bxd,axd,impuls);
pxd=pxd/max(pxd);


