% Finite-difference solution of the 2-D cochlear model
% Neely (1981, JASA 69, 1386-1393)


function fdm81()
global r1 r2 r4 srd y m;
L=3.5;
%% Admittanz
figure(1);clf;
subplot(211);hold on;axis([0 L -130 -30]);
title('Admittanz');
ylabel('Amplitude in dB');
subplot(212);hold on;axis([0 L -1 1]);
xlabel('Abstand Steigbügel in cm');
ylabel('Phase in \pi');

%% Druck
figure(2);clf;
subplot(211);hold on;axis([0 L 70 170]);
title('Druck');
ylabel('Amplitude in dB');
subplot(212);hold on;axis([0 L -9 1]);
xlabel('Abstand Steigbügel in cm');
ylabel('Phase in \pi');

%% Auslenkung
figure(3);clf;
subplot(211);hold on;axis([0 L -70 30]);
title('Auslenkung BM');
ylabel('Amplitude in dB');
subplot(212);hold on;axis([0 L -9 1]);
xlabel('Abstand Steigbügel in cm');
ylabel('Phase in \pi');
drawnow;

%% Parameterwerte
m=8; n=246; xl=3.5; yh=0.1; rho=1;
k0=1e9; r0=200; m0=0.15; k1=-2;

%% Array definieren
a = zeros(m,m,n);
p = zeros(m,n);

%% Variablen definieren
dx=xl/(n-1);
dy=yh/(m-1);
r1=(dx/dy)^2;
r2=r1+r1;
r4=2+r2;

%% Schleife mit 5 Frequenzen
for fi=1:5,
   w=2*pi*200*(2^fi);
   s=i*w;
   srd=i*4*rho*w*dy*r1;
   p0p=2*rho*w^2;
   
   % Berechnung der Admittanz für alle x
   x=transpose(linspace(0,xl,n));
   y=(k0*exp(k1*x)/s + r0 + m0*s).^(-1);
   
   % Lösung der Blockmatrizen
   %  Durchgang von oben nach unten
  	p(:,1)=-2*dx*p0p;   
   a(:,:,1)=mxadd(1);
	b=inv(reshape(a(:,:,1),m,m));
	a(:,:,1)=reshape(-2*b,1,m,m);
  	p(:,1)=b*p(:,1);   
   for k=2:n-1,
      a(:,:,k)=a(:,:,k-1)+mxadd(k);
		b=inv(reshape(a(:,:,k),m,m));
		a(:,:,k)=reshape(-b,1,m,m);
  		p(:,k)=b*p(:,k-1);   
   end
   %  Durchgang von unten nach oben
   for k=n-2:-1:1,
      b=reshape(a(:,:,k),m,m);
		p(:,k)=p(:,k)-b*p(:,k+1);
   end
   d=y.*p(1,:).'./s;
   
   % Datensatz plotten
   list{fi} = [num2str(w/(2*pi)) ' Hz'];
   legend(list,'Location','southoutside', 'orientation', 'horizontal');
   figure(1);
   subplot(211);plot(x,dbmag(y));
   subplot(212);plot(x,phase(y));
   legend(list,'Location','southoutside', 'orientation', 'horizontal');
   figure(2);
   subplot(211);plot(x,dbmag(p(1,:)));
   subplot(212);plot(x,phase(p(1,:)));
   legend(list,'Location','southoutside', 'orientation', 'horizontal');
   figure(3);
   subplot(211);plot(x,dbmag(d));
   subplot(212);plot(x,phase(d));
   legend(list,'Location','southoutside', 'orientation', 'horizontal');
   drawnow;

end

%% Algorithmus
function aa=mxadd(k)
global r1 r2 r4 srd y m;
aa=zeros(m,m,1);
aa(1,1,1)=r4+y(k)*srd;
aa(1,2,1)=-r2;
for j=2:m-1,
   aa(j,j-1,1)=-r1;
   aa(j,j,1)=r4;
   aa(j,j+1,1)=-r1;
end
aa(m,m-1,1)=-r2;
aa(m,m,1)=r4;

%% Transformation von Phase und Amplitude in dB
function y=dbmag(x)
y=20*real(log10(max(1e-9,x)));

function y=phase(x)
y=unwrap(imag(log(max(1e-9,x))))/pi;