% Time Doppler plot
function [logTD,time,doppl]=tdoppl(nc,txpol,rangebin,xpix,ypix)

if nargin<3, xpix=250; end
if nargin<4, ypix=80; end

% load I and Q data
[I,Q]=ipixcdf(nc,txpol,rangebin);
N=size(I,1);

% adjust fft window and averaging to match image size
wdw=max(128,2^ceil(log(4*N/xpix)/log(2)));
yavg=max(1,ceil(wdw/ypix));
timeStep=wdw/4;
M=floor((N-wdw)/timeStep);

% short time fourier transforms
hw=hamming(wdw);
TD=zeros(floor(wdw/yavg),M);
y=zeros(yavg,floor(wdw/yavg)); leny=prod(size(y));
for m=1:M,
  i=(m-1)*timeStep;
  x=I(1+i:wdw+i)+j*Q(1+i:wdw+i);
  x=hw.*(I(1+i:wdw+i)+j*Q(1+i:wdw+i));
  fftx=abs(fftshift(fft(x)));
  y(:)=fftx(1:leny);
  TD(:,m)=mean(y,1)';
end

% time and normalized frequency
PRF=nc{'PRF'}(1);              % Pulse Repetition Frequency [Hz]
time=[0 M]/PRF;
freq=[-0.5 0.5]*PRF;

% convert to doppler velocity
doppl=freq*3e8/(2*nc{'RF_frequency'}(1)*1e9); 

% enhance log-plot with noise floor
mn=mean(TD,2);
[mn,indx]=sort(mn);
noise=TD(indx(1:2),:);
noiseFloor=median(noise(:));
TD(find(TD<noiseFloor))=noiseFloor;
logTD=log(TD);
mn=min(min(logTD)); mx=max(max(logTD));
logTD=(logTD-mn)*63/(mx-mn);

% display image
if nargout<1,
  time=(wdw/2+(0:M-1)*timeStep)/PRF;
  freq=((0:wdw/yavg-1)/(wdw/yavg)-0.5)*PRF;
  doppl=freq*3e8/(2*nc{'RF_frequency'}(1)*1e9); 
  image(time,doppl,logTD); set(gca,'ydir','normal');
end
