% thltdr.m - reflectance and diameter from EMAV THL files
%
function thltdr
global fmax xmax T D fdw srf mbw wd lcr pac pep smz
fmax=20;    % maximum frequency (kHz)
xmax=12;    % maximum tube length (cm)
D=0.8;      % tube diameter (cm)
T=25;       % air temperature (deg C)
fdw=1;      % frequency-domain window ?
srf=4;      % sampling-rate factor
mbw=18;     % maximum bandwidth of frequency-domain window (kHz)
pac=0;      % plot autocorrelation ?
pep=0;      % plot epsilon(x) & B(x) ?
lcr=1;      % compute cavity length from reflectance ?
smz=1;      % smooth load impedance ?
wd=0;       % write data ?
%
fnlst=get_fn('*.THL');
if (isempty(fnlst))
   fprintf('No THL files.\n');
   return
end
[set,fnlst]=select_set(fnlst,['P','C']);
plot_set(fnlst,set); % fd refl, td refl, & diameter plots
return

function [set,fnlst]=select_set(fnlst,chr)
cp=6;
nc=length(chr);
nf=length(fnlst);
kk=zeros(nf:nc);
for k=1:nc
   for j=1:nf
      kk(j,k)=(fnlst{j}(cp)==chr(k));
   end
end
ns=sum(kk);
nn=sum(ns>0);
if (nn<1)
   set=0;
elseif (nn==1)
   set=find(ns>0);
else
   for k=1:nc
      fprintf('%5d %s\n',k,chr(k));
   end
   set=input(sprintf('\nWhich file set (1-%d)? [0] ',nc));
   if (isempty(set)) set=0; end
end
if (set>0)
   fnlst=fnlst(logical(kk(:,set)));
end
return

function filelist=get_fn(fp)
fn = dir(fp);
nfn=length(fn);
if (nfn < 1)
   filelist = [];
   return; 
end
for k=1:nfn
    filelist(k) = {fn(k).name};
end
return

% get cavity pressures from THL file
function [f,pl,zl,sr]=load_thl(fn,ch)
if (isoctave)
   load(fn);         % fetch data
else
   load(fn,'-MAT');  % fetch data
end
sr=rate;
if (ch == 1)
	pl=[pl1(:);mean(pl1)];
	zl=[zl1(:);mean(zl1)];
else
	pl=[pl2(:);mean(pl2)];
	zl=[zl2(:);mean(zl2)];
end
n=length(zl)-1;
f=(0:n)'*df;
return

function plot_set(fnlst,set)
global pac pep
nfig=3;
if (pep) nfig=4; end
for fig=1:nfig
   figure(fig);clf
end
for i=1:length(fnlst)
   if (pac) figure(9);clf; end
   plot_refl(char(fnlst(i)),set,i);
   drawnow
end
figure(1)
subplot(2,1,1);hold off
subplot(2,1,2);hold off
if (nfig>1) figure(2); hold off; end
if (nfig>2) figure(3); hold off; end
if (nfig>3) figure(4); hold off; end
return

function plot_refl(fn,set,seq)
global fmax xmax T fdw srf mbw wd lcr pac pep smz
c=plt_color(seq);
[f,pl1,zl1,rate]=load_thl(fn,1);
[f,pl2,zl2,rate]=load_thl(fn,2);
if (smz)
   zo1=surge(zl1,z_tube,fdw);
   zo2=surge(zl2,z_tube,fdw);
   zl1=zsmo(zl1,zo1,rate);
   zl2=zsmo(zl2,zo2,rate);
end
f=f/1000;
fmax=min([fmax rate/2000]);
flim=[1/fmax fmax];
if (fdw>0)
   if (mbw<(rate/2000))
      fdw=mbw*2000/rate;
   else
      fdw=1;
   end
end
%
n=length(zl1)-1;
lc1=cavity_length(pl1,rate)*10;
lc2=cavity_length(pl2,rate)*10;
[zo1,zi1,zx1]=surge(zl1,z_tube,fdw);
[zo2,zi2,zx2]=surge(zl2,z_tube,fdw);
R1=reflect(zl1,zx1);
R2=reflect(zl2,zx2);
fprintf('%s: ',fn);
fprintf('len= %4.1f %4.1f; ',lc1,lc2);
fprintf('srg= %5.1f %5.1f; ',zo1,zo2);
fprintf('nrt= %5.1f %5.1f\n',zi1,zi2);
if (pac) plot_autocor(pl1,pl2,lc1,lc2,rate,fn); end
%
figure(1)
dblim=[-8.5 0.5];
phlim=[-6 1];
gdlim=[-0.2 1];
db1=20*log10(abs(R1)+eps);
db2=20*log10(abs(R2)+eps);
ph1=unwrap(angle(R1))/(2*pi);
ph2=unwrap(angle(R2))/(2*pi);
gd1=delay(R1,f);
gd2=delay(R2,f);
subplot(2,1,1);
semilogx(f,db1,c)
axis([flim dblim]);
ylabel('magnitude (dB)')
title('frequency-domain reflectance');
hold on
subplot(2,1,2);
semilogx(f,gd1,c)
axis([flim gdlim]);
ylabel('delay (ms)')
xlabel('frequency (kHz)')
hold on
%
subplot(2,1,1);
semilogx(f,db2,c)
hold on
subplot(2,1,2);
semilogx(f,gd2,c)
hold on
drawnow
%
figure(2)
lim=[-0.1 1.1 -4 16];
r1=ffs(fd_window(R1,srf,fdw));
r2=ffs(fd_window(R2,srf,fdw));
r1=fftshift(r1);
r2=fftshift(r2);
n=length(r1)/2;
rk=rate*srf/1000;
t=((-n):(n-1))/rk;
plot(t,r1*rk,c)
axis(lim)
title('time-domain reflectance');
xlabel('time (ms)')
ylabel('(kHz)')
hold on
plot(t,r2*rk,c)
hold on
drawnow
if (lcr)
   lc1=cavlen_tdr(t,r1);
   lc2=cavlen_tdr(t,r1);
end
r1=fftshift(r1);
r2=fftshift(r2);
%
figure(3)
ur=rate*srf;
[A1,B1,x]=refl_inverse(r1,zo1,ur,xmax,T);
[A2,B2,x]=refl_inverse(r2,zo2,ur,xmax,T);
D1=20*sqrt(A1/pi);
D2=20*sqrt(A2/pi);
E1=[diff(log(max(A1,eps)))./diff(x)/2 0];
E2=[diff(log(max(A2,eps)))./diff(x)/2 0];
x=x*10;
lim=[-5 75 -0.5 16.5];
plot(x,D1,c)
hold on
plot(lc1,interp1(x,D1,lc1),strcat(c,'o'))
axis(lim)
title('inverse solution');
xlabel('distance (mm)')
ylabel('diameter (mm)')
plot(x,D2,c)
hold on
plot(lc2,interp1(x,D2,lc2),strcat(c,'o'))
%
if (pep)
   figure(4)
   cc=sprintf('%s--',c);
   plot(x,B1,c)
   title('logarithmic gradient');
   xlabel('distance (mm)')
   axis([-5 75 -4 1])
   hold on
   plot(x,E1,cc)
   plot(x,B2,c)
   plot(x,E2,cc)
   hold on
end
%
drawnow
if (wd)
   write_data(sprintf('rf1_%d%d.txt',set,seq),[f db1 gd1]);
   write_data(sprintf('rf2_%d%d.txt',set,seq),[f db2 gd2]);
   write_data(sprintf('rt1_%d%d.txt',set,seq),[t(:) r1(:)]);
   write_data(sprintf('rt2_%d%d.txt',set,seq),[t(:) r2(:)]);
   write_data(sprintf('dx1_%d%d.txt',set,seq),[x(:) D1(:)]);
   write_data(sprintf('dx2_%d%d.txt',set,seq),[x(:) D2(:)]);
end
return

% inverse solution for radius from reflectance
function [A,B,x]=refl_inverse(r,zo,rate,xmax,T);
rho = 1.1769e-3 * (1 - 0.00335 * (T - 26.85));
c   = 3.4723e4  * (1 + 0.00166 * (T - 26.85));
dt=1/rate;
dx=c*dt;
Ao=rho*c/zo;
m=ceil(xmax/dx);
if (length(r)<(2*m))
   m=floor(length(r)/2);
   xmax=m*dx;
   fprintf('reduced xmax=%5.2f\n',xmax);
end
x=(0:(m-1))*dx;
epmx=1/dx;      % maximum epsilon
etmx=1/dx;      % maximum eta
pfmx=1.01;      % maximum power flow
pfmn=-0.01;     % minimum power flow
%
n=2*m;
pm=r(1:n);
pp=zeros(n,1);
pp(1)=1;
pm(1)=0;
A=zeros(1,m);
B=zeros(1,m);
A(1)=Ao;
for k=1:(m-1)
   p=pp+pm;                             % pressure
   q=cumsum(p)*dt;                      % integrated pressure
   u=pp-pm+c*B(k)*q;                    % velocity / Yo
   a=A(k)/Ao;                           % relative area
   pf=p(k)*conj(u(k))*a;                % power flow
   if ((pf>pfmx)|(pf<pfmn)) break; end  % check #1
   j1=k;
   j2=k+1;
   bk=B(k)/2;
   v1=bk*(p(j1)+u(j1))-pm(j1)/dx;
   v2=bk*(p(j2)+u(j2))-pm(j2)/dx;
   dd=u(j1)*q(j2)-q(j1)*u(j2);
   if (dd==0) break; end                % check #2
   ep=(q(j2)*v1-q(j1)*v2)/dd;
   et=(u(j1)*v2-u(j2)*v1)/dd;
   if ((ep>epmx)) break; end            % check #3
   if ((et>etmx)&(bk<-1)) break; end    % check #4
   npp=zeros(n,1);
   npm=zeros(n,1);
   for j=(k+1):(n-k)
      jp=j+1;
      jm=j-1;
      vp=bk*p(jp)-((ep-bk)*u(jp)+et*q(jp));
      vm=bk*p(jm)+((ep-bk)*u(jm)+et*q(jm));
      npp(j)=pp(jm)-vm*dx;
      npm(j)=pm(jp)-vp*dx;
   end
   pp=npp;
   pm=npm;
   A(k+1)=A(k)*exp(2*ep*dx);
   B(k+1)=B(k)+2*et*dt;
end
if (0)
   fprintf('x=%5.2f A=%6.3g B=%6.3g ',x(k),A(k),B(k));
   fprintf('ep=%6.3g et=%6.3g pf=%6.3g\n',ep,et,pf);
end;
return

function plot_autocor(pc1,pc2,lc1,lc2,rate,fn)
global T
c=3.4723e4*(1+0.00166*(T-26.85));
a1=ffs(abs(pc1).^2);
a2=ffs(abs(pc2).^2);
a1=a1/a1(1);
a2=a2/a2(1);
n=length(a1);
t=(0:(n-1))*1000/rate;
figure(9)
plot(t,a1,t,a2)
hold on
xlim([0 0.5])
title(fn)
xlabel('time (ms)')
ylabel('pressure autocorrelation')
t1=200*lc1/c;
t2=200*lc2/c;
x1=[t1 t1];
x2=[t2 t2];
y1=[min(a1) max(a1)];
y2=[min(a2) max(a2)];
plot(x1,y1,':',x2,y2,':')
return

% cavity length load pressure autocorrelation
function [lk,td]=cavity_length(pc,rate)
global T
lmn=1;      % minimum length (cm)
lmx=12;     % maximum length (cm)
c=3.4723e4*(1+0.00166*(T-26.85));
%
imn=round(lmn*rate*2/c)+1;
imx=round(lmx*rate*2/c)+1;
a=ffs(abs(pc).^2);
while (a(imn+1)<a(imn)) imn=imn+1; end % increase imn
[amx,m]=max(a(imn:imx));
if ((amx/a(1))>0.5)
   m=m+imn-1;
else
   m=imn+1;
	while (a(m+1)>a(m)) m=m+1; end
end
d=(a(m-1)-a(m+1))/(a(m-1)-2*a(m)+a(m+1))/2;
td=(m+d-1)/rate;
lk=td*c/2;
return

% cavity length based on TDR peak
function lc = cavlen_tdr(t,r)
global T
c = 3.4723e4 * (1 + 0.00166 * (T - 26.85));
[mx,ii]=max(r);
for k=ii:length(r); if (r(k)<(mx/2)) break; end; end
lc=c*t(k)/200;
return

% smooth load impedance by restricting TDR via LR
function zl=zsmo(zl,zo,rate)
fr=[0.020 20];   % frequency range (kHz)
ti=[-0.5   2];   % time interval (ms)
nf=length(zl);
nt=2*(nf-1);
rk=rate/1000;
fr=round(fr/rk*nt);
ti=round(ti*rk);
f=(fr(1):fr(2))';
t=(ti(1):ti(2))/nt;
ii=f+1;
A=exp(-i*2*pi*f*t);
AA=[real(A);imag(A)];
rl=(zl(ii)-zo)./(zl(ii)+zo);
rt=AA\[real(rl);imag(rl)];
rf=zeros(size(zl));
rf(ii)=A*rt;
zl=zo*(1+rf)./(1-rf);
return

function R=reflect(zl,zx)
R=(zl-zx)./(zl+conj(zx));
return

% group delay
function gd=delay(R,f)
[n,m] = size(R);
ph = unwrap(angle(R))/(2*pi);
gd = zeros([n,m]);
for k=1:m
   gd(:,k) = -cdif(ph(:,k))./cdif(f(:));
end
return

% centered difference
function dx=cdif(x)
n=length(x);
dx=zeros(size(x));
dx(1)=x(2)-x(1);
dx(2:(n-1))=(x(3:n)-x(1:(n-2)))/2;
dx(n)=x(n)-x(n-1);
return

% characteristic impedance
function zo = z_tube
global T D
c = 3.4723e4 * (1 + 0.00166 * (T - 26.85));
rho = 1.1769e-3 * (1 - 0.00335 * (T - 26.85));
zo = (rho * c) / (pi * (D / 2)^2);
return

% frequency-domain window
% usage: H = fd_window(H,srf,fdw)
% H   - transfer function
% srf - sampling-rate factor
% fdw - bandwidth / Nyquist_frequency
function H=fd_window(H,srf,fdw)
if (fdw<=0) return; end
n=length(H)-1;
p=pi*(0:n)'/n/fdw;
a=0.16;                 % Blackman window
w=(1-a+cos(p)+a*cos(2*p))/2;
w(p>pi)=0;
Z=zeros((srf-1)*n,1);   % zero padding
H=[H.*w;Z];
return

% surge impedance 
% usage: [zo,zi,zx]=surge(zl,zo,fdw)
% zl  - load impedance
% zo  - tube impedance
% fdw - bandwidth / Nyquist_frequency
function [zo,zi,zx]=surge(zl,zo,fdw)
nf=length(zl);
w=2*pi*(0:(nf-1))'/(nf-1);
s=i*w;
v0=fd_window(ones(size(zl)),1,fdw);
v1=fd_window(sin(w),1,fdw);
v0=v0/sum(v0);
v1=v1/sum(v1);
z1=-sum(imag(zl).*v1);
zn=z1/2;
zi=abs(sum(imag(zl).*w)/sum(w.*w))/4;
zx=zo+zn+zi*s;
for k=1:40        % iterate
   rl=(zl-zx)./(zl+conj(zx));
   zo=zo*(1+sum(real(rl).*v0));
   zi=zi*(1+sum(imag(rl).*v1));
   zx=zo+zn+zi.*s;
end
return

function c=plt_color(seq)
colors=['b' 'g' 'r' 'c' 'm' 'k'];
c=colors(mod(seq-1,6)+1);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fast Fourier synthesize real signal
function h=ffs(H)
m=length(H);
n=2*(m-1);
H(1,:)=real(H(1,:));
H(m,:)=real(H(m,:));
H((m+1):n,:)=conj(H((m-1):-1:2,:));
h=real(ifft(H));
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function o = isoctave
o=1;eval('OCTAVE_VERSION;','o=0;');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function write_data(fn,data)
[nr,nc] = size(data);
fp=fopen(fn,'wt');
fprintf(fp,'; %s\n', fn);
for i=1:nr
    for j=1:nc;
        fprintf(fp,' %14.5g',data(i,j));
    end
    fprintf(fp,'\n');
end
fclose(fp);