Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/AuditoryToolbox/SeneffEarSetup.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [SeneffPreemphasis, SeneffFilterBank, SeneffForward, SeneffBackward] ...
wolffd@0 2 = SeneffEarSetup(fs)
wolffd@0 3
wolffd@0 4 % This m-function is based on data from the following paper:
wolffd@0 5 % Benjamin D. Bryant and John D. Gowdy, "Simulation of Stages
wolffd@0 6 % I and II of Seneff's Auditory Model (SAM) Using Matlab", and
wolffd@0 7 % published in the Proceedings of the 1993 Matlab User's Group
wolffd@0 8 % Conference.
wolffd@0 9 % Thanks to Benjamin Bryant for supplying us with his filter
wolffd@0 10 % coefficients and the initial organization of this implementation.
wolffd@0 11
wolffd@0 12 % (c) 1998 Interval Research Corporation
wolffd@0 13
wolffd@0 14 % Set the following variable to a non-zero value to see a summary
wolffd@0 15 % of the filter bank's behaviour.
wolffd@0 16 plotTests = 0;
wolffd@0 17
wolffd@0 18 % The following values were taken from Figure 2 of Bryant's paper.
wolffd@0 19 PreemphasisRTheta = [0.86 3.1148863;0.99 0; 0.5 0; 0.95 3.14159];
wolffd@0 20
wolffd@0 21 % The following values were taken from Table 1 of Bryant's paper.
wolffd@0 22 % They represent the cascade zeros (R-z and Theta-z), and the
wolffd@0 23 % second order poles (radius and theta) and zeros (radius and theta/2).
wolffd@0 24 %
wolffd@0 25 % R-z Theta-z Radius Theta R-z2
wolffd@0 26 FilterBankRTheta = [
wolffd@0 27 0 3.14159 0.740055 2.633909 0.8
wolffd@0 28 0.86 2.997077 0.753637 2.178169 0.8
wolffd@0 29 0.86 2.879267 0.775569 1.856744 0.8
wolffd@0 30 0.86 2.761458 0.798336 1.617919 0.8
wolffd@0 31 0.86 2.643648 0.819169 1.433496 0.8
wolffd@0 32 0.86 2.525839 0.837158 1.286795 0.8
wolffd@0 33 0.8 2.964876 0.852598 1.167321 0.8
wolffd@0 34 0.86 2.408029 0.865429 1.068141 0.8
wolffd@0 35 0.86 2.29022 0.876208 0.984489 0.8
wolffd@0 36 0.86 2.17241 0.885329 0.912985 0.8
wolffd@0 37 0.86 2.054601 0.893116 0.851162 0.8
wolffd@0 38 0.86 1.936791 0.899823 0.797179 0.8
wolffd@0 39 0.8 2.788161 0.906118 0.749633 0.8
wolffd@0 40 0.86 1.818981 0.911236 0.70744 0.8
wolffd@0 41 0.86 1.701172 0.915747 0.669742 0.8
wolffd@0 42 0.86 1.583362 0.919753 0.635858 0.8
wolffd@0 43 0.86 1.465552 0.923335 0.605237 0.8
wolffd@0 44 0.86 1.347743 0.926565 0.57743 0.8
wolffd@0 45 0.8 2.611447 0.929914 0.552065 0.8
wolffd@0 46 0.86 1.229933 0.932576 0.528834 0.8
wolffd@0 47 0.86 1.112123 0.944589 0.487783 0.75
wolffd@0 48 0.86 0.994314 0.957206 0.452645 0.660714
wolffd@0 49 0.86 0.876504 0.956548 0.42223 0.672143
wolffd@0 50 0.86 0.758694 0.956653 0.395644 0.682143
wolffd@0 51 0.8 2.434732 0.956518 0.372208 0.690966
wolffd@0 52 0.86 0.640885 0.956676 0.351393 0.69881
wolffd@0 53 0.86 0.523075 0.956741 0.316044 0.712143
wolffd@0 54 0.8 2.258018 0.956481 0.287157 0.723052
wolffd@0 55 0.8 2.081304 0.956445 0.263108 0.732143
wolffd@0 56 0.8 1.904589 0.956481 0.242776 0.739835
wolffd@0 57 0.86 0.405265 0.958259 0.217558 0.749384
wolffd@0 58 0.8 1.727875 0.963083 0.197086 0.757143
wolffd@0 59 0.8 1.55116 0.969757 0.175115 0.769048
wolffd@0 60 0.8 1.374446 0.97003 0.153697 0.780662
wolffd@0 61 0.8 1.197732 0.970382 0.134026 0.791337
wolffd@0 62 0.8 1.021017 0.970721 0.118819 0.799596
wolffd@0 63 0.8 1.5 0.970985 0.106711 0.8
wolffd@0 64 0.8 1.2 0.971222 0.096843 0.8
wolffd@0 65 0.8 1 0.97144 0.088645 0.8
wolffd@0 66 0.8 0.9 0.971645 0.081727 0.8];
wolffd@0 67
wolffd@0 68 % Let's plot the cascade zero locations and the locations of the
wolffd@0 69 % pole and zeros in the resonator.
wolffd@0 70 if plotTests
wolffd@0 71 clf;
wolffd@0 72 subplot(3,3,1);
wolffd@0 73 plot(FilterBankRTheta(:,1).*exp(i*FilterBankRTheta(:,2)))
wolffd@0 74 axis([-1 1 0 1])
wolffd@0 75 title('Cascade Zero Locations')
wolffd@0 76
wolffd@0 77 subplot(3,3,2);
wolffd@0 78 plot([FilterBankRTheta(:,3).*exp(i*FilterBankRTheta(:,4)) ...
wolffd@0 79 FilterBankRTheta(:,5).*exp(i*FilterBankRTheta(:,4)/2)],'+')
wolffd@0 80 title('Resonator Pole/Zero')
wolffd@0 81 drawnow;
wolffd@0 82 end
wolffd@0 83
wolffd@0 84 % Convert r-theta form, first into a list of roots, then a polynomial
wolffd@0 85 roots=exp(i*PreemphasisRTheta(:,2)).*PreemphasisRTheta(:,1);
wolffd@0 86 SeneffPreemphasis=real(poly([roots;conj(roots)]));
wolffd@0 87
wolffd@0 88 % Plot the preemphasis filter response, if desired
wolffd@0 89 if plotTests
wolffd@0 90 subplot(3,3,3);
wolffd@0 91 freqScale=(0:255)/256*8000;
wolffd@0 92 freqresp = FreqResp(SeneffPreemphasis,[1], freqScale, 16000);
wolffd@0 93 semilogx(freqScale,freqresp)
wolffd@0 94 title('Preemphasis Response');
wolffd@0 95 axis([100 10000 -60 20])
wolffd@0 96 drawnow;
wolffd@0 97 end
wolffd@0 98
wolffd@0 99 % Now figure out the second order sections that make up the main
wolffd@0 100 % filter bank cascade. We put the zeros into the numerator (b's)
wolffd@0 101 % and there are no poles. Just to keep things simpler, we adjust
wolffd@0 102 % the gain of each filter to keep it unity gain at DC.
wolffd@0 103 [channels,width] = size(FilterBankRTheta);
wolffd@0 104 roots=exp(i*FilterBankRTheta(:,2)).*FilterBankRTheta(:,1);
wolffd@0 105 SeneffFilterBank = zeros(channels,5);
wolffd@0 106 for j=1:channels
wolffd@0 107 SeneffFilterBank(j,1:3) = poly([roots(j) conj(roots(j))]);
wolffd@0 108 SeneffFilterBank(j,1:3) = SeneffFilterBank(j,1:3)/sum(SeneffFilterBank(j,1:3));
wolffd@0 109 end
wolffd@0 110
wolffd@0 111 % Plot the cascade zero responses, if desired.
wolffd@0 112 if plotTests
wolffd@0 113 subplot(3,3,4);
wolffd@0 114 y=soscascade([1 zeros(1,511)],SeneffFilterBank);
wolffd@0 115 freqresp=20*log10(abs(fft(y(1:5:40,:)')));
wolffd@0 116 freqScale=(0:511)/512*16000;
wolffd@0 117 semilogx(freqScale(1:256),freqresp(1:256,:))
wolffd@0 118 axis([100 10000 -150 0]);
wolffd@0 119 title('Cascade Response');
wolffd@0 120 drawnow;
wolffd@0 121 end
wolffd@0 122
wolffd@0 123 % Now figure out the resonating filters. Each of these resonators
wolffd@0 124 % is a double pole-zero pair.
wolffd@0 125 zlocs = FilterBankRTheta(:,5).*exp(i*FilterBankRTheta(:,4)/2);
wolffd@0 126 plocs = FilterBankRTheta(:,3).*exp(i*FilterBankRTheta(:,4));
wolffd@0 127 SeneffForward = zeros(5,channels);
wolffd@0 128 SeneffBackward = zeros(5,channels);
wolffd@0 129
wolffd@0 130 for j=1:channels
wolffd@0 131 SeneffForward(:,j) = real(poly([zlocs(j) conj(zlocs(j)) ...
wolffd@0 132 zlocs(j) conj(zlocs(j))]))';
wolffd@0 133 SeneffBackward(:,j) = real(poly([plocs(j) conj(plocs(j)) ...
wolffd@0 134 plocs(j) conj(plocs(j))]))';
wolffd@0 135 end
wolffd@0 136
wolffd@0 137 % Now plot the frequency response of just the resonating filters.
wolffd@0 138 % These are all bandpass filters.
wolffd@0 139 if plotTests
wolffd@0 140 subplot(3,3,5);
wolffd@0 141 impulse = [1 zeros(1,255)];
wolffd@0 142 y=zeros(256,channels);
wolffd@0 143 for j=1:40
wolffd@0 144 y(:,j) = filter(SeneffForward(:,j),SeneffBackward(:,j),impulse)';
wolffd@0 145 end
wolffd@0 146 freqresp=20*log10(abs(fft(y(:,1:5:40))));
wolffd@0 147 freqScale=(0:255)/256*16000;
wolffd@0 148 semilogx(freqScale(1:128),freqresp(1:128,:))
wolffd@0 149 axis([100 10000 -30 40]);
wolffd@0 150 title('Resonators Response')
wolffd@0 151 drawnow;
wolffd@0 152 end
wolffd@0 153
wolffd@0 154 % The plot below shows the overall response of the preemphasis filters
wolffd@0 155 % along with the just-designed cascade of zeros.
wolffd@0 156 if plotTests
wolffd@0 157 subplot(3,3,6);
wolffd@0 158 impulse = [1 zeros(1,511)];
wolffd@0 159 y=soscascade(filter(SeneffPreemphasis, [1], impulse), ...
wolffd@0 160 SeneffFilterBank);
wolffd@0 161 freqresp=20*log10(abs(fft(y(1:5:40,:)')));
wolffd@0 162 freqScale=(0:511)/512*16000;
wolffd@0 163 semilogx(freqScale(1:256),freqresp(1:256,:))
wolffd@0 164 axis([100 10000 -100 25]);
wolffd@0 165 title('Preemphasis+Cascade');
wolffd@0 166 drawnow;
wolffd@0 167 end
wolffd@0 168
wolffd@0 169 % Now we need to normalize the gain of each channel. We run an impulse
wolffd@0 170 % through the preemphasis filter, and then through the cascade of zeros.
wolffd@0 171 % Finally, we run it through each of the resonator filters.
wolffd@0 172 impulse = [1 zeros(1,255)];
wolffd@0 173 y=soscascade(filter(SeneffPreemphasis, [1], impulse), ...
wolffd@0 174 SeneffFilterBank);
wolffd@0 175 for j=1:channels
wolffd@0 176 y(j,:) = filter(SeneffForward(:,j),SeneffBackward(:,j),y(j,:));
wolffd@0 177 end
wolffd@0 178
wolffd@0 179 % Now we have impulse responses from each filter. We can find the FFT
wolffd@0 180 % and then find the gain peak. We divide each forward polynomial by the
wolffd@0 181 % maximum gain (to normalize) and then multiply by the desired low
wolffd@0 182 % frequency roll-off. The Bryant paper says that the last 24 channels
wolffd@0 183 % should be cut at 6dB per octave and that this occurs at 1600 Hz, but
wolffd@0 184 % it looks to me like the gain change happens at 3200 Hz.
wolffd@0 185 freqresp=abs(fft(y'));
wolffd@0 186 gain = ones(1,channels)./max(freqresp);
wolffd@0 187 cfs = FilterBankRTheta(:,4)/pi*fs/2;
wolffd@0 188 rolloff = min(cfs/1600,1);
wolffd@0 189
wolffd@0 190 for j=1:channels
wolffd@0 191 SeneffForward(:,j)=SeneffForward(:,j)*gain(j)*rolloff(j);
wolffd@0 192 end
wolffd@0 193
wolffd@0 194 % All Done. The figure below should match Figure 3 of Bryant's paper.
wolffd@0 195 if plotTests
wolffd@0 196 subplot(3,3,8);
wolffd@0 197 impulse = [1 zeros(1,511)];
wolffd@0 198 y=soscascade(filter(SeneffPreemphasis, [1], impulse), ...
wolffd@0 199 SeneffFilterBank);
wolffd@0 200 for j=1:channels
wolffd@0 201 y(j,:) = filter(SeneffForward(:,j),SeneffBackward(:,j),y(j,:));
wolffd@0 202 end
wolffd@0 203
wolffd@0 204 freqresp=20*log10(abs(fft(y(1:5:40,:)')));
wolffd@0 205 freqScale=(0:511)/512*16000;
wolffd@0 206 plot(freqScale(1:256),freqresp(1:256,:))
wolffd@0 207 axis([100 10000 -120 0]);
wolffd@0 208 title('Magnitude Response vs. Linear Frequency');
wolffd@0 209 drawnow;
wolffd@0 210 end
wolffd@0 211
wolffd@0 212
wolffd@0 213 function mag=FreqResp(b,a,f,fs)
wolffd@0 214 cf = exp(i*2*pi*f/fs);
wolffd@0 215 num = 0;
wolffd@0 216 for i=1:length(b)
wolffd@0 217 num = num + b(end-i+1)*cf.^i;
wolffd@0 218 end
wolffd@0 219
wolffd@0 220 denom = 0;
wolffd@0 221 for i=1:length(a)
wolffd@0 222 denom = denom + a(end-i+1)*cf.^i;
wolffd@0 223 end
wolffd@0 224 mag = 20*log10(abs(num./denom));