function codesonar(fs,sb2,cpb,cod,range)
% codesonar(fs,sb2,cpb,cod,range)
% Long-range sonar using binary codes.
% Necessary hardware: Speaker connected to the left audio output (up to 20 kHz),
% Microphone at the left audio input (up to 20 kHz)
% fs     : sample rate                         (Default: 48000)
% sb2    : samples per bit = 2^sb2             (Default: 7, so 128 samples per bit)
% cpb    : number of carrier cycles per bit    (Default: 2)
% cod    : binary code                         (Default: Barker code of 13 bits)
% range  : maximum range in meters             (Default: 300)

% K. v. d. Heide,  e-mail:  v.d.heide@on-line.de
%
% Parameter-Defaults
% ==================
if nargin<5, range = 300;   end           % maximum range in meters
if nargin<4,   cod = barkercode(13)>0; end% binary code
if nargin<3,   cpb = 2;     end           % number of carrier cycles per bit
if nargin<2,   sb2 = 7;     end           % samples per bit = 2^sb2
if nargin<1,    fs = 48000; end           % sampling rate in Hz

% Constants
% =========
cs    = 340;                              % speed of sound in m/s
tcy   = 2*(fs/cs*(range+2)+length(cod)*2^sb2)/fs;
if     tcy<0.1, tc = 0.1;
elseif tcy<0.5, tc = 0.05;
else            tc = 0.0;
end
tcs   = round(fs*tc);                     % samples to receive within tc
spb   = 2^sb2;                            % samples per bit
f0    = cpb/spb;                          % carrier frequency relative to fs
ptn   = length(cod)*spb;                  % length of pulsetrain in samples
nc    = round(2*(range+2)/cs*fs) + ptn;   % #samples to receive
car   = exp(-1i*2*pi*f0*(1:nc)');         % analytic carrier
cds   = nfach(2*cod-1,spb)';              % code in samples
cdr   = reverse(2*cod-1)';                % reverse code for correlation
pulse = cds.*imag(car(1:size(cds,1),1));  % pulse signal
txsig = [pulse; zeros(nc+tcs,1)];         % tx signal in complete cycle
q     = min(0.4,2*f0);
gh    = 2*ceil(1/q);                      % necessary degree of Hilbert filter
bh    = firpm(gh,[q,1-q],[1 1],'Hilbert');% design of Hilbertfilter
relax = 0.01;                             % relaxationfactor for mean echo
mps   = cs/fs/2*2^sb2;                    % meter per sample decimated
minre = ptn/fs*cs;                        % minimum distance
n     = sb2;
fnts  = 14;                               % FontSize
cycl  = nc + ptn + tcs;                   % cycle in samples
bfr   = ceil(1*fs/cycl);                  % output buffer in complete cycles

% Graphic
% =======
mi = nc; 
fb = 1:mi;
fh = fullscreen('Code Sonar','w',3);
lh = line((0:mi-1)*mps,zeros(1,mi));      % echo line
lm = line((0:mi-1)*mps,zeros(1,mi),'Color','r');  % mean echo line
ax = gca;
set(ax,'XLim',[minre range],'YLim',[0 1.25])
xlabel('distance [m]','FontSize',fnts)
ylabel('echo amplitude','FontSize',fnts)
title('Code-SONAR','FontSize',fnts)
cdstr = '01';
text(minre+0.1,1.20,['sample rate:       ' num2str(fs)],   'FontSize',fnts,'FontName','Courier')
text(minre+0.1,1.15,['carrier frequency: ' num2str(f0*fs)],'FontSize',fnts,'FontName','Courier')
text(minre+0.1,1.10,['pulse:             ' cdstr(cod+1)],  'FontSize',fnts,'FontName','Courier')
text(minre+0.1,1.05,['range:             ' num2str(minre) ' - ' num2str(range) ' m'], ...
                     'FontSize',fnts,'FontName','Courier')

% 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 = nc+ptn+tcs;
ai.TimeOut       = 1.0;
ao = analogoutput('winsound');
ch = addchannel(ao,1);
ao.BufferingMode = 'auto';
ao.SampleRate    = fs;

% Start
% =====
cnt = 0;
x   = cell(n+1);                          % signal structure
xs  = complex(zeros(1,n));                % last sample
for k=1:n+1, x{k} = complex(zeros(ptn*2^(k-1),1)); end
for k=1:bfr
   putdata(ao,txsig)
end
flushdata(ai)
start([ao ai])                            % start sampling

% Period
% ======
while 1
  try
      
      % sampling
      putdata(ao,txsig);                  % output tx signal
      rxs = getdata(ai,cycl);             % wait for cycle input samples
      rxs = rxs((ptn+1):end-tcs);
      
      % Hilbert filter
      rxa  = rxs + 1i*delay(filter(bh,1,rxs),-gh/2); % Hilbert filter
      
      % sfift to baseband
      x{1} = car.*rxa;                    % spectral rotation into baseband
      
      % multirate lowpass filter
      for k=1:n                           % simple multirate lowpass filter
         x{k+1} = 0.5*x{k}(1:2:end-1) + 0.25*([xs(k); x{k}(2:2:end-2)] + x{k}(2:2:end));
         xs(k)  = x{k}(end);              % save last sample for next loop cycle
      end
      
      % correlation
      yd = abs(filter(cdr,1,x{n+1}));     % correlation
      
      % graphic
      if cnt==0
         [~,ix] = max(yd);
         ly = length(yd);
         y  = [yd(ix:end); zeros(ix-1,1)];
         ym = relax*y;
      else
         y  = [yd(ix:end); zeros(ix-1,1)];
         ym = relax*y + (1-relax)*ym;     % mean echo
      end
      yy = max(yd(ix:ly));
      yym = max(ym);
      set(lh,'XData',(1:ly)*mps,'YData',y/yy)   % update echo
      set(lm,'XData',(1:ly)*mps,'YData',ym/yym) % update mean echo
      drawnow
      cnt = cnt + 1;                      % time count in seconds
   catch
      return
   end
end
