% fdm86.m - compute frequency-domain cochlear model 
% from Neely & Kim (1986, JASA 79, 1472-1480).
function fdm86()
%
% misc constants
%
g=1;
b=0.4;
gamma=1;
rho=1;
nf=4;
%
% spatial dimensions
%
L=2.5;
H=0.1;
W=0.1;
N=251;
%
% arrays
%
x=transpose(linspace(0,L,N));
o=ones(N-1,1);
P=zeros(N,nf);
Yb=zeros(N,nf);
Db=zeros(N,nf);
Dh=zeros(N,nf);
%
% BM stifness, damping, and mass
%
k1=1.1e9*exp(-4*x);
c1=20+1500*exp(-2*x);
m1=0.003;
k2=7e6*exp(-4.4*x);
c2=10*exp(-2.2*x);
m2=0.0005;
k3=1e7*exp(-4*x);
c3=2*exp(-0.8*x);
k4=6.15e8*exp(-4*x);
c4=1040*exp(-2*x);
%
% Middle-ear stifness, damping, and mass
%
km=2.1e5;
cm=400;
mm=0.045;
As=0.01;
Am=0.35;
Gm=0.5;
Pe=2.848e-4;
%
% loop over frequencies
%
f=100;
dx=L/(N-1);
Ap=b/g;
for J=1:nf,
   f=f*4;
   s=2i*pi*f;
   srd=2*s*rho*dx;
   % impedance
   z1=k1/s+c1+m1*s;
   z2=k2/s+c2+m2*s;
   z3=k3/s+c3;
   z4=k4/s+c4;
   Zb=z1+z2.*(z3-gamma*z4)./(z2+z3);
   Zp=Zb/Ap;
   hc=z2./(z2+z3);
   % middle ear
   Zm=km/s+cm+mm*s;
   ame=-srd*As/Zm;
   pme=ame*(Am/(Gm*As))*Pe;
   % solve for pressure
   Q=zeros(N,1);
   Q(1)=pme;
   a=(s*rho*dx*dx/H)*(Zp.^(-1));
   A=diag(2+a,0)+diag(-o,1)+diag(-o,-1);
   A(1,1)=1+ame;
   p=A\Q;
   P(:,J)=p/Pe;
   Yb(:,J)=Zb.^(-1);
   Db(:,J)=((p./Zb)./s)*1e7;
   Dh(:,J)=hc.*Db(:,J);
end
%
% plot figures
%
xx=x/L;
tpi=2*pi;
% impedance
figure(1);clf;
subplot(211);plot(xx,real(20*log10(Yb)));axis([0 1 -120 -20]);
title('impedance');
subplot(212);plot(xx,unwrap(imag(log(Yb)))/tpi);axis([0 1 -1 1]);
xlabel('x/L');
% pressure
figure(2);clf;
subplot(211);plot(xx,real(20*log10(P)));axis([0 1 -40 60]);
title('pressure');
subplot(212);plot(xx,unwrap(imag(log(P)))/tpi);axis([0 1 -5 1]);
xlabel('x/L');
% BM displacement
figure(3);clf;
subplot(211);plot(xx,real(20*log10(Db)));axis([0 1 -100 0]);
title('BM displacement');
subplot(212);plot(xx,unwrap(imag(log(Db)))/tpi);axis([0 1 -5 1]);
xlabel('x/L');
% HC displacement
figure(4);clf;
subplot(211);plot(xx,real(20*log10(Dh)));axis([0 1 -100 0]);
title('HC displacement');
subplot(212);plot(xx,unwrap(imag(log(Dh)))/tpi);axis([0 1 -5 1]);
xlabel('x/L');