function w = chirpsonar(fs,fu,fo,tc,rw,mde,buffer,TM)
% w = chirpsonar(fs,fu,fo,tc,rw,mde,buffer,TM)
% Sonar with linear chirp (sawtooth) between edge frequencies fu and fo (in Hz).
% Necessary hardware: Speaker connected to the left audio output (up to 20 kHz),
% Microphone at the left audio input (up to 20 kHz)
% fs     : sampling rate in Hz                      (Default: 96000)
% fu     : lower edge frequency of the chirp in Hz  (Default: 12000)
% fo     : upper edge frequency of the chirp in Hz  (Default: 20000)
% tc     : cycle time in seconds                    (Default: 1/3 s)
% rw     : range in meters (upper limit of the plot (Default: 6)
% mde=0  : mode of the display is sattic, otherwise dynamic (sensitive to change, Default: 1)
% buffer : output buffer in cycles                  (Default: 3)
% TM     : total time of the experiment in seconds  (Default: inf, abort by deleting figure)

% K. v. d. Heide,  e-mail:  v.d.heide@on-line.de
%
% Parameter-Defaults
% ==================
if nargin<8, TM = inf;   end              % Duration of experiment in seconds
if nargin<7, buffer = 3; end              % sze of output buffer in cycles
if nargin<6, mde = 1;    end              % display mode
if nargin<5, rw = 6;     end              % sonar range in meters
if nargin<4, tc = 1/3;   end              % cycle time in seconds
if nargin<3, fo = 20000; end              % upper edge frequency of the chirp
if nargin<2, fu = 12000; end              % lower edge frequency of the chirp
if nargin<1, fs = 96000; end              % sampling rate in Hz
cs = 340;                                 % spedd of sound in m/s
if rw > 40*tc, error('rw > 40*tc not allowed'); end
if fo<=fu,     error('fO <= fu');               end
if fo>fs/2,    error('fo > Nyquist frequency'); end

% Constants
% =========
g     = 70;                               % order of Hilbertfilter
q     = min(0.2,fu/(fs/2));               % passband of Hilbertfilter
bh    = firpm(g,[q,1-q],[1 1],'Hilbert'); % Hilbertfilter
fnts  = 14;                               % FontSize
n     = round(fs*tc);                     % blocklength in Samples
tc    = n/fs;                             % actual cycle time 
dd    = 2*(fo-fu)*rw/cs/n;                % maximum frequency shift (rel. Nyquist)
[b,a] = butter(4,dd);                     % IIR-filter design
hgn   = hamming(n)';                      % Hamming window
relax = 0.1;                              % relaxation factor for mean echo
ns0   = 3;

% Graphic
% =======
tr = 2*rw/cs;                             % Laufzeit bei Reichweite rw
ft = (fo-fu)/tc;                          % df/dt in Hz/s
df = ft*tr;                               % frequency shift in Hz
mi = 2*round(n*df/fs) + 500;              % index into FFT at maximum range
nn = g/2 + 2;                             % delay in samples
fb = [n-nn+1:n  1:mi-nn];                 % indexvektor into FFT
m_pro_bin = cs*tc*fs/(2*n*(fo-fu));       % mapping factor bin ==> meters
af = (1:mi).^4;                           % amplitude factor over distance
fh = fullscreen('Sonar','w',3);
lh = line((0:mi-1)*m_pro_bin,zeros(1,mi));% plot echo (abs(fft))
ls = line((0:mi-1)*m_pro_bin,zeros(1,mi),'Color','r','LineWidth',2);% plot mean echo
ax = gca;
set(ax,'XLim',[0 rw])
xlabel('distance [m]','FontSize',fnts)
ylabel('echo signal','FontSize',fnts)
title('Chirp-SONAR: Echo Amplitude over Distnce','FontSize',fnts)

% Delete all running DAQs
% =======================
daqopen = daqfind;
for k=1:length(daqopen)
   stop(daqopen(k));
   delete(daqopen(k));
end

% Define Audio Channels
% =====================
ai = analoginput('winsound');
addchannel(ai,1);
ai.SampleRate    = fs;
ai.TriggerRepeat = inf;
ai.TriggerType   = 'immediate';
ai.SamplesPerTrigger = n;
ai.TimeOut       = 2*tc;
ao = analogoutput('winsound');
ch = addchannel(ao,1:2);
set(ao,'SampleRate',fs)
ao.TimeOut       = 2*tc;
flushdata(ai)                             % delete old input data

% Start
% =====
ph    = 2*pi*(fu+(fo-fu)/tc/2*(0:n-1)/n*tc).*(0:n-1)/n*tc;  % chirp phase
txa   = exp(1i*ph);                       % complex chirp for shift to base band
txsig = real(txa)';                       % 
uu    = 2;
txsig(1:uu) = txsig(1:uu).*(1-cos(pi/(uu+1)*(1:uu)'));
txsig(end-(0:uu-1)) = txsig(end-(0:uu-1)).*(1-cos(pi/(uu+1)*(1:uu)'));
txdat = [txsig  zeros(n,1)];              % chirp signal on left channel, zeros on right
for k=1:buffer+1
   putdata(ao,txdat);                     % output chirp data into output buffer
end
start([ai ao])                            % start input and output
rxsig = getdata(ai,n)';                   % wait for n input samples
T     = 0;

% Period
% ======
yx = 0;
ym = zeros(1,mi);
while T<TM
   try
      rxsig = getdata(ai,n)';             % wait for n input samples
      putdata(ao,txdat);                  % output chirp data 
      rxa = rxsig(1:n) - 1i*delay(filter(bh,1,rxsig(1:n)),-g/2);
      v   = txa.*rxa;                     % spectral rotation into baseband
      w   = filter(b,a,v);                % lowpass filter 
      sp  = fft(w.*hgn);                  % FFT
      yd  = af.*abs(sp(fb));              % amplitude scaling by distance^4
      ym  = relax*yd + (1-relax)*ym;      % mean echo
      [~,ix] = max(ym);
      if mde
         y = yd - ym;
      else
         y = yd;
      end
      ys  = y(ix+ns0:end);
      yy  = max(ys);
      if yy>yx,     yx = 2*yy; end 
      if yy<0.2*yx, yx = 2*yy; end 
      set(lh,'XData',(1:length(ys))*m_pro_bin,'YData',ys)             % update echo
      set(ls,'XData',(1:length(ys))*m_pro_bin,'YData',ym(ix+ns0:end)) % update mean echo
      set(ax,'YLim',[0 yx])
      drawnow
      T  = T + tc;                        % time count in seconds
   catch
      return
   end
end
