%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Signal Processing Simulation of Melbourne Nucleus Cochlear Implant %
%                                                                    %
% Katherine J. Herrick and Jane Ueda                                 %
%                                                                    %
% Based on " A New Portable Sound Processor for the University of    % 
% Melbourne/Nuclues Limited Multielectrode Cochlear Implant"         %
% by Hugh J. McDermontt, Colette M. McKay, and Andrew E. Vandali     %
% (J. Acoust. Soc. Am. 91(6), June 1992)                             %
%                                                                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;				% clear all variables

channel=16;			% number of channels
load dilmbin			% load sound data
srate=8192;			% sound sampling rate
N=length(wavedata);		% length of sound record
dt=1/srate;			% sampling interval
df=1./(dt.*N);			% frequency interval
fmax=df.*N./2;			% maximum frequency
t=(1: N)*dt-dt;			% time array
f=(1: N)*df-df;			% frequency array
f1=(0:511)*3400/512;		% frequency array for sampled spectrum
Tmax=N*dt;			% time length of sound record
d=.5*srate;			% frequency scalar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Time Waveform of Original Sound %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
subplot(2,1,1),plot(t, wavedata);
title('Original Time Waveform of "Dilemma"');
xlabel('Time (sec)');		% label x axis
ylabel('Amplitude');		% label y axis
axis([0 Tmax -3e4 3e4]);	% set axis limits

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Spectrogram of Sound %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,1,2),specgram(wavedata,[],8192,[],[]);
title('Spectrogram of Original Sound');


%%%%%%%%%%%%%%%%%%%%%%%
% Play original sound %
%%%%%%%%%%%%%%%%%%%%%%%

%sound(wavedata, 8192);		% First time
pause(.1)
%sound(wavedata, 8192);		% Second time

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define bandwidths for channels %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

w(1)=200;

for i=1: channel,
   w(i+1)=(w(i)+200);
end   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Filter sound into 16 channels using  %
% 8th Order Chebyshev Bandpass Filters %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[b1, a1]=cheby2(3,30,[w(1)/d w(2)/d]);		% First channel
[b2, a2]=cheby2(3,30,[w(2)/d w(3)/d]);		% Second channel
[b3, a3]=cheby2(3,30,[w(3)/d w(4)/d]);		% Third channel
[b4, a4]=cheby2(3,30,[w(4)/d w(5)/d]);		% Fourth channel
[b5, a5]=cheby2(3,30,[w(5)/d w(6)/d]);		% Fifth channel
[b6, a6]=cheby2(3,30,[w(6)/d w(7)/d]);		% Sixth channel
[b7, a7]=cheby2(3,30,[w(7)/d w(8)/d]);		% Seventh channel
[b8, a8]=cheby2(3,30,[w(8)/d w(9)/d]);		% Eighth channel
[b9, a9]=cheby2(3,30,[w(9)/d w(10)/d]);		% Ninth channel
[b10, a10]=cheby2(3,30,[w(10)/d w(11)/d]);		% Tenth channel
[b11, a11]=cheby2(3,30,[w(11)/d w(12)/d]);		% Eleventh channel
[b12, a12]=cheby2(3,30,[w(12)/d w(13)/d]);		% Twelfth channel
[b13, a13]=cheby2(3,30,[w(13)/d w(14)/d]);		% Thriteenth channel
[b14, a14]=cheby2(3,30,[w(14)/d w(15)/d]);		% Fourteenth channel
[b15, a15]=cheby2(3,30,[w(15)/d w(16)/d]);		% Fifteenth channel
[b16, a16]=cheby2(3,30,[w(16)/d w(17)/d]);		% Sixteenth channel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Z transform Channel Signals %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h1 = freqz(b1, a1, 512, srate);		% 1st channel
h2 = freqz(b2, a2, 512, srate);		% 2nd channel
h3 = freqz(b3, a3, 512, srate);		% 3rd channel
h4 = freqz(b4, a4, 512, srate);		% 4th channel
h5 = freqz(b5, a5, 512, srate);		% 5th channel
h6 = freqz(b6, a6, 512, srate);		% 6th channel
h7 = freqz(b7, a7, 512, srate);		% 7th channel
h8 = freqz(b8, a8, 512, srate);		% 8th channel
h9 = freqz(b4, a4, 512, srate);		% 9th channel
h10 = freqz(b10, a10, 512, srate);		% 10th channel
h11 = freqz(b11, a11, 512, srate);		% 11th channel
h12 = freqz(b12, a12, 512, srate);		% 12th channel
h13 = freqz(b13, a13, 512, srate);		% 13th channel
h14 = freqz(b14, a14, 512, srate);		% 14th channel
h15 = freqz(b15, a15, 512, srate);		% 15th channel
h16 = freqz(b16, a16, 512, srate);		% 16th channel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Frequency Spectrum of First and Sixteenth Channel % 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)
subplot(2,1,1), plot(f1, abs(h1));	%Plot first channel
title('16 Chebyshev bandpass filters allocate the time-waveform data into bins.');

subplot(2,1,2), plot(f1, abs(h16));	%Plot sixteenth channel
title('Displayed are the first(200-400 Hz) and the sixteenth(3200-3400 Hz)');
pause(2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fourth Order Lowpass Chebyshev Filter %
% Cutoff Frequency= 200 Hz              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[blow, alow]=cheby2(4,30,200/d);

%%%%%%%%%%%%%%%%%%
%Impulse
x=zeros(32,1);
x(1)=1;
x=x';
out(1,:)=filter(b1, a1, wavedata);
out(2,:)=filter(b2, a2, wavedata);
out(3,:)=filter(b3, a3, wavedata);
out(4,:)=filter(b4, a4, wavedata);
out(5,:)=filter(b5, a5, wavedata);
out(6,:)=filter(b6, a6, wavedata);
out(7,:)=filter(b7, a7, wavedata);
out(8,:)=filter(b8, a8, wavedata);
out(9,:)=filter(b9, a9, wavedata);
out(10,:)=filter(b10, a10, wavedata);
out(11,:)=filter(b11, a11, wavedata);
out(12,:)=filter(b12, a12, wavedata);
out(13,:)=filter(b13, a13, wavedata);
out(14,:)=filter(b14, a14, wavedata);
out(15,:)=filter(b15, a15, wavedata);
out(16,:)=filter(b16, a16, wavedata);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Time Waveform of Channel Signal %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)
subplot(2,1,1), plot(t, out(1,:));
title('Filtered Time Waveform');
ylabel('1st Channel: 200-400 Hz');
xlabel('Time (s)');
axis([0 2.2 -8000 8000]);
subplot(2,1,2), plot(t, out(16,:));
ylabel('16th Channel: 3200-3400 Hz');
xlabel('Time (s)');
axis([0 2.2 -8000 8000]);
pause(2);


%%%%%%%%%%%%%%%%%%%%%%%%
% Signal Rectification %
%%%%%%%%%%%%%%%%%%%%%%%%

  for i=1:channel, 
    out(i,:)= abs(out(i,:));
    low(i,:)=filter(blow, alow, out(i,:));
  end
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Rectified Time Waveform of Channel Signal %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(4)
subplot(2,1,1), plot(t, abs(out(1,:)));
title('Rectified Filtered Time Waveform');
ylabel('1st Channel: 200-400 Hz');
xlabel('Time (s)');
axis([0 2.2 -8000 8000]);

subplot(2,1,2), plot(t, abs(out(16,:)));
ylabel('16th Channel: 3200-3400 Hz');
xlabel('Time (s)');
axis([0 2.2 -8000 8000]);
pause(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Envelope of Rectified Time Waveform of Channel Signal %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(5)
subplot(2,1,1), plot(t, low(1,:));
title('Envelope of Rectified Filtered Time Waveform');
ylabel('1st 200-400 Hz filter');
xlabel('Time (s)');
axis([0 2.2 -8000 8000]);
subplot(2,1,2), plot(t, low(16,:));
ylabel('16th 3200-3400 Hz filter');
xlabel('Time (s)');
axis([0 2.2 -8000 8000]);
pause(2);
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Impulse Response of Channel Filters %
% Used for synthesis of sound         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 vec(1, 1:32)=filter(b1, a1, x);
 vec(2, 1:32)=filter(b2, a2, x);
 vec(3, 1:32)=filter(b3, a3, x);
 vec(4, 1:32)=filter(b4, a4, x);
 vec(5, 1:32)=filter(b5, a5, x);
 vec(6, 1:32)=filter(b6, a6, x);
 vec(7, 1:32)=filter(b7, a7, x);
 vec(8, 1:32)=filter(b8, a8, x);
 vec(9, 1:32)=filter(b9, a9, x);
 vec(10, 1:32)=filter(b10, a10, x);
 vec(11, 1:32)=filter(b11, a11, x);
 vec(12, 1:32)=filter(b12, a12, x);
 vec(13, 1:32)=filter(b13, a13, x);
 vec(14, 1:32)=filter(b14, a14, x);
 vec(15, 1:32)=filter(b15, a15, x);
 vec(16, 1:32)=filter(b16, a16, x);

%%%%%%%%%%%%%%%%%%%%%%%%%%
% Maximum Peak Selection % 
%%%%%%%%%%%%%%%%%%%%%%%%%%
 for n=1:523				% For total number of windows
   sample=[low(1:16, n*32-16)];		% Sample channel signal for window 
   sample=sample';			% Transpose vector of samples
   [s, k]=sort(sample);			% Sort channel strengths 

   % Scale and sum six strongest channels
   output(n*32-31: n*32)=s(11:16)*vec(k(11:16),:); 

 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Time Waveform of Synthesized Sound %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(6);
subplot(2,1,1),plot(t(1:16736), output);
title('Synthesized Time Waveform of "Dilemma"');
xlabel('Time (sec)');
ylabel('Amplitude');
axis([0 Tmax -200 200]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plot Spectrogram of Synthesized sound %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 subplot(2,1,2),specgram(output,[],8192,[],[]);
title('Spectrogram of Synthesized Sound'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% playback synthesized sound %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pause(.1);
%sound(output, 8192);
pause(0.1);
%sound(output, 8192);

%%%%%%%%%%%%%%%%%%%%%
% End of Simulation %
%%%%%%%%%%%%%%%%%%%%%