wolffd@0: function [SeneffPreemphasis, SeneffFilterBank, SeneffForward, SeneffBackward] ...
wolffd@0: 			= SeneffEarSetup(fs)
wolffd@0: 			
wolffd@0: % This m-function is based on data from the following paper:
wolffd@0: %	Benjamin D. Bryant and John D. Gowdy, "Simulation of Stages
wolffd@0: %	I and II of Seneff's Auditory Model (SAM) Using Matlab", and 
wolffd@0: %	published in the Proceedings of the 1993 Matlab User's Group
wolffd@0: %	Conference.
wolffd@0: % Thanks to Benjamin Bryant for supplying us with his filter 
wolffd@0: % coefficients and the initial organization of this implementation.
wolffd@0: 
wolffd@0: % (c) 1998 Interval Research Corporation  
wolffd@0: 
wolffd@0: % Set the following variable to a non-zero value to see a summary
wolffd@0: % of the filter bank's behaviour.
wolffd@0: plotTests = 0;
wolffd@0: 
wolffd@0: % The following values were taken from Figure 2 of Bryant's paper.
wolffd@0: PreemphasisRTheta = [0.86 3.1148863;0.99 0; 0.5 0; 0.95 3.14159];
wolffd@0: 
wolffd@0: % The following values were taken from Table 1 of Bryant's paper.
wolffd@0: % They represent the cascade zeros (R-z and Theta-z), and the 
wolffd@0: % second order poles (radius and theta) and zeros (radius and theta/2).
wolffd@0: %
wolffd@0: %    R-z       Theta-z       Radius       Theta       R-z2
wolffd@0: FilterBankRTheta = [
wolffd@0:      0         3.14159       0.740055     2.633909    0.8
wolffd@0:      0.86      2.997077      0.753637     2.178169    0.8
wolffd@0:      0.86      2.879267      0.775569     1.856744    0.8
wolffd@0:      0.86      2.761458      0.798336     1.617919    0.8
wolffd@0:      0.86      2.643648      0.819169     1.433496    0.8
wolffd@0:      0.86      2.525839      0.837158     1.286795    0.8
wolffd@0:      0.8       2.964876      0.852598     1.167321    0.8
wolffd@0:      0.86      2.408029      0.865429     1.068141    0.8
wolffd@0:      0.86      2.29022       0.876208     0.984489    0.8
wolffd@0:      0.86      2.17241       0.885329     0.912985    0.8
wolffd@0:      0.86      2.054601      0.893116     0.851162    0.8
wolffd@0:      0.86      1.936791      0.899823     0.797179    0.8
wolffd@0:      0.8       2.788161      0.906118     0.749633    0.8
wolffd@0:      0.86      1.818981      0.911236     0.70744     0.8
wolffd@0:      0.86      1.701172      0.915747     0.669742    0.8
wolffd@0:      0.86      1.583362      0.919753     0.635858    0.8
wolffd@0:      0.86      1.465552      0.923335     0.605237    0.8
wolffd@0:      0.86      1.347743      0.926565     0.57743     0.8
wolffd@0:      0.8       2.611447      0.929914     0.552065    0.8
wolffd@0:      0.86      1.229933      0.932576     0.528834    0.8
wolffd@0:      0.86      1.112123      0.944589     0.487783    0.75
wolffd@0:      0.86      0.994314      0.957206     0.452645    0.660714
wolffd@0:      0.86      0.876504      0.956548     0.42223     0.672143
wolffd@0:      0.86      0.758694      0.956653     0.395644    0.682143
wolffd@0:      0.8       2.434732      0.956518     0.372208    0.690966
wolffd@0:      0.86      0.640885      0.956676     0.351393    0.69881
wolffd@0:      0.86      0.523075      0.956741     0.316044    0.712143
wolffd@0:      0.8       2.258018      0.956481     0.287157    0.723052
wolffd@0:      0.8       2.081304      0.956445     0.263108    0.732143
wolffd@0:      0.8       1.904589      0.956481     0.242776    0.739835
wolffd@0:      0.86      0.405265      0.958259     0.217558    0.749384
wolffd@0:      0.8       1.727875      0.963083     0.197086    0.757143
wolffd@0:      0.8       1.55116       0.969757     0.175115    0.769048
wolffd@0:      0.8       1.374446      0.97003      0.153697    0.780662
wolffd@0:      0.8       1.197732      0.970382     0.134026    0.791337
wolffd@0:      0.8       1.021017      0.970721     0.118819    0.799596
wolffd@0:      0.8       1.5           0.970985     0.106711    0.8
wolffd@0:      0.8       1.2           0.971222     0.096843    0.8
wolffd@0:      0.8       1             0.97144      0.088645    0.8
wolffd@0:      0.8       0.9           0.971645     0.081727    0.8];
wolffd@0: 
wolffd@0: % Let's plot the cascade zero locations and the locations of the
wolffd@0: % pole and zeros in the resonator.
wolffd@0: if plotTests
wolffd@0: 	clf;
wolffd@0: 	subplot(3,3,1);
wolffd@0: 	plot(FilterBankRTheta(:,1).*exp(i*FilterBankRTheta(:,2)))
wolffd@0: 	axis([-1 1 0 1])
wolffd@0: 	title('Cascade Zero Locations')
wolffd@0: 
wolffd@0: 	subplot(3,3,2);
wolffd@0: 	plot([FilterBankRTheta(:,3).*exp(i*FilterBankRTheta(:,4)) ...
wolffd@0: 		  FilterBankRTheta(:,5).*exp(i*FilterBankRTheta(:,4)/2)],'+')
wolffd@0: 	title('Resonator Pole/Zero')
wolffd@0: 	drawnow;
wolffd@0: end
wolffd@0: 
wolffd@0: % Convert r-theta form, first into a list of roots, then a polynomial
wolffd@0: roots=exp(i*PreemphasisRTheta(:,2)).*PreemphasisRTheta(:,1);
wolffd@0: SeneffPreemphasis=real(poly([roots;conj(roots)]));
wolffd@0: 
wolffd@0: % Plot the preemphasis filter response, if desired
wolffd@0: if plotTests
wolffd@0: 	subplot(3,3,3);
wolffd@0: 	freqScale=(0:255)/256*8000;
wolffd@0: 	freqresp = FreqResp(SeneffPreemphasis,[1], freqScale, 16000);
wolffd@0: 	semilogx(freqScale,freqresp)
wolffd@0: 	title('Preemphasis Response');
wolffd@0: 	axis([100 10000 -60 20])
wolffd@0: 	drawnow;
wolffd@0: end
wolffd@0: 
wolffd@0: % Now figure out the second order sections that make up the main
wolffd@0: % filter bank cascade.  We put the zeros into the numerator (b's)
wolffd@0: % and there are no poles.  Just to keep things simpler, we adjust
wolffd@0: % the gain of each filter to keep it unity gain at DC.
wolffd@0: [channels,width] = size(FilterBankRTheta);
wolffd@0: roots=exp(i*FilterBankRTheta(:,2)).*FilterBankRTheta(:,1);
wolffd@0: SeneffFilterBank = zeros(channels,5);
wolffd@0: for j=1:channels
wolffd@0: 	SeneffFilterBank(j,1:3) = poly([roots(j) conj(roots(j))]);
wolffd@0: 	SeneffFilterBank(j,1:3) = SeneffFilterBank(j,1:3)/sum(SeneffFilterBank(j,1:3));
wolffd@0: end
wolffd@0: 
wolffd@0: % Plot the cascade zero responses, if desired.
wolffd@0: if plotTests
wolffd@0: 	subplot(3,3,4);
wolffd@0: 	y=soscascade([1 zeros(1,511)],SeneffFilterBank);
wolffd@0: 	freqresp=20*log10(abs(fft(y(1:5:40,:)')));
wolffd@0: 	freqScale=(0:511)/512*16000;
wolffd@0: 	semilogx(freqScale(1:256),freqresp(1:256,:))
wolffd@0: 	axis([100 10000 -150 0]);
wolffd@0: 	title('Cascade Response');
wolffd@0: 	drawnow;
wolffd@0: end
wolffd@0: 
wolffd@0: % Now figure out the resonating filters.  Each of these resonators
wolffd@0: % is a double pole-zero pair.
wolffd@0: zlocs = FilterBankRTheta(:,5).*exp(i*FilterBankRTheta(:,4)/2);
wolffd@0: plocs = FilterBankRTheta(:,3).*exp(i*FilterBankRTheta(:,4));
wolffd@0: SeneffForward = zeros(5,channels);
wolffd@0: SeneffBackward = zeros(5,channels);
wolffd@0: 
wolffd@0: for j=1:channels
wolffd@0: 	SeneffForward(:,j) = real(poly([zlocs(j) conj(zlocs(j)) ...
wolffd@0: 										zlocs(j) conj(zlocs(j))]))';
wolffd@0: 	SeneffBackward(:,j) = real(poly([plocs(j) conj(plocs(j)) ...
wolffd@0: 										plocs(j) conj(plocs(j))]))';
wolffd@0: end
wolffd@0: 
wolffd@0: % Now plot the frequency response of just the resonating filters.
wolffd@0: % These are all bandpass filters.
wolffd@0: if plotTests
wolffd@0: 	subplot(3,3,5);
wolffd@0: 	impulse = [1 zeros(1,255)];
wolffd@0: 	y=zeros(256,channels);
wolffd@0: 	for j=1:40
wolffd@0: 		y(:,j) = filter(SeneffForward(:,j),SeneffBackward(:,j),impulse)';
wolffd@0: 	end
wolffd@0: 	freqresp=20*log10(abs(fft(y(:,1:5:40))));
wolffd@0: 	freqScale=(0:255)/256*16000;
wolffd@0: 	semilogx(freqScale(1:128),freqresp(1:128,:))
wolffd@0: 	axis([100 10000 -30 40]);
wolffd@0: 	title('Resonators Response')
wolffd@0: 	drawnow;
wolffd@0: end
wolffd@0: 
wolffd@0: % The plot below shows the overall response of the preemphasis filters
wolffd@0: % along with the just-designed cascade of zeros.
wolffd@0: if plotTests
wolffd@0: 	subplot(3,3,6);
wolffd@0: 	impulse = [1 zeros(1,511)];
wolffd@0: 	y=soscascade(filter(SeneffPreemphasis, [1], impulse), ...
wolffd@0: 				SeneffFilterBank);
wolffd@0: 	freqresp=20*log10(abs(fft(y(1:5:40,:)')));
wolffd@0: 	freqScale=(0:511)/512*16000;
wolffd@0: 	semilogx(freqScale(1:256),freqresp(1:256,:))
wolffd@0: 	axis([100 10000 -100 25]);
wolffd@0: 	title('Preemphasis+Cascade');
wolffd@0: 	drawnow;
wolffd@0: end
wolffd@0: 
wolffd@0: % Now we need to normalize the gain of each channel.  We run an impulse
wolffd@0: % through the preemphasis filter, and then through the cascade of zeros.
wolffd@0: % Finally, we run it through each of the resonator filters.
wolffd@0: impulse = [1 zeros(1,255)];
wolffd@0: y=soscascade(filter(SeneffPreemphasis, [1], impulse), ...
wolffd@0: 			SeneffFilterBank);
wolffd@0: for j=1:channels
wolffd@0: 	y(j,:) = filter(SeneffForward(:,j),SeneffBackward(:,j),y(j,:));
wolffd@0: end
wolffd@0: 
wolffd@0: % Now we have impulse responses from each filter.   We can find the FFT
wolffd@0: % and then find the gain peak.  We divide each forward polynomial by the
wolffd@0: % maximum gain (to normalize) and then multiply by the desired low
wolffd@0: % frequency roll-off.  The Bryant paper says that the last 24 channels
wolffd@0: % should be cut at 6dB per octave and that this occurs at 1600 Hz, but 
wolffd@0: % it looks to me like the gain change happens at 3200 Hz.
wolffd@0: freqresp=abs(fft(y'));
wolffd@0: gain = ones(1,channels)./max(freqresp);
wolffd@0: cfs = FilterBankRTheta(:,4)/pi*fs/2;
wolffd@0: rolloff = min(cfs/1600,1);
wolffd@0: 	
wolffd@0: for j=1:channels
wolffd@0: 	SeneffForward(:,j)=SeneffForward(:,j)*gain(j)*rolloff(j);
wolffd@0: end
wolffd@0: 	
wolffd@0: % All Done. The figure below should match Figure 3 of Bryant's paper.
wolffd@0: if plotTests
wolffd@0: 	subplot(3,3,8);
wolffd@0: 	impulse = [1 zeros(1,511)];
wolffd@0: 	y=soscascade(filter(SeneffPreemphasis, [1], impulse), ...
wolffd@0: 				SeneffFilterBank);
wolffd@0: 	for j=1:channels
wolffd@0: 		y(j,:) = filter(SeneffForward(:,j),SeneffBackward(:,j),y(j,:));
wolffd@0: 	end
wolffd@0: 
wolffd@0: 	freqresp=20*log10(abs(fft(y(1:5:40,:)')));
wolffd@0: 	freqScale=(0:511)/512*16000;
wolffd@0: 	plot(freqScale(1:256),freqresp(1:256,:))
wolffd@0: 	axis([100 10000 -120 0]);
wolffd@0: 	title('Magnitude Response vs. Linear Frequency');
wolffd@0: 	drawnow;
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: function mag=FreqResp(b,a,f,fs)
wolffd@0: cf = exp(i*2*pi*f/fs);
wolffd@0: num = 0;
wolffd@0: for i=1:length(b)
wolffd@0: 	num = num + b(end-i+1)*cf.^i;
wolffd@0: end
wolffd@0: 
wolffd@0: denom = 0;
wolffd@0: for i=1:length(a)
wolffd@0: 	denom = denom + a(end-i+1)*cf.^i;
wolffd@0: end
wolffd@0: mag = 20*log10(abs(num./denom));