Mercurial > hg > map

comparison userProgramsASRforDummies/cEssexAid.m @ 38:c2204b18f4a2 tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
End nov big change
author Ray Meddis <rmeddis@essex.ac.uk>
date Mon, 28 Nov 2011 13:34:28 +0000
parents
children
comparison
equal deleted inserted replaced
37:771a643d5c29 38:c2204b18f4a2
1 classdef cEssexAid
2 %ESSEXAID_WRAPCLASS Wrapper for the EssexAid - Nick Clark July 2011
3 % This class wraps up the EssexAid algorithm function that processes
4 % each block of samples. This wrapper closely emulates the GUI used
5 % in the lab and runs stimuli through the exact same algorithm used
6 % in the lab. It even includes a helper function to generate C code
7 % from the algorithm for use in a real-time framework.
8
9
10 %% *********************************************************
11 % properties _ _
12 % | | (_)
13 % _ __ _ __ ___ _ __ ___ _ __| |_ _ ___ ___
14 % | '_ \| '__/ _ \| '_ \ / _ \ '__| __| |/ _ \/ __|
15 % | |_) | | | (_) | |_) | __/ | | |_| | __/\__ \
16 % | .__/|_| \___/| .__/ \___|_| \__|_|\___||___/
17 % | | | |
18 % |_| |_|
19 %************************************************************
20
21 %% **********************************************************
22 % Public properties - can be set by user
23 %************************************************************
24 properties(Access = public)
25 sr = 48e3;
26 numSamples = 1024; %MAX=6912, LAB_USE=48
27 stimulusUSER
28
29 %------------------------------------------------------------------
30 % Params for audiometric freqs 250, 500, 1000, 2000, 4000, 8000 Hz
31 %------------------------------------------------------------------
32 audiometry_dB= [ 0; 0; 0; 0; 0; 0]; %Pure tone threshold in dB SPL
33 mainGain_dB = [ 0; 0; 0; 0; 0; 0]; %Gain applied at audiometric frequencies
34 TC_dBHL = [40; 40; 40; 40; 40; 40]; %Compression thresholds (in dB HL from 2nd filt)
35 TM_dBHL = [10; 10; 10; 10; 10; 10]; %MOC thresholds (in dB OUTPUT from 2nd filt)
36 DRNLc = [ 0.2; 0.2; 0.2; 0.2; 0.2; 0.2]; %Compression exponent at audiometric frequencies
37
38 %------------------------------------------------------------------
39 % Dynamic compression properties
40 %------------------------------------------------------------------
41 ARtau = 60e-3; %decay time constant
42 ARthreshold_dB = 85; %dB SPL (input signal level) =>200 to disable
43 MOCtau = 450e-3; %Time constant in Seconds
44 MOCfactor = 0.5; %dB attenuation applied to the input per dB exceeding output threshold
45
46 %------------------------------------------------------------------
47 % Band filtering properties
48 %------------------------------------------------------------------
49 bwOct = 1/2; %1/1, 1/2, 1/3, 1/4, 1/5
50 filterOrder = 2 %BUTTER=2, GTF=3
51 useGTF = false; %If false, revert to butterworth
52 end
53
54 %% **********************************************************
55 % Read only properties that are not dependent
56 %************************************************************
57 properties(SetAccess = private)
58 MOCrecord
59 end
60
61 %% **********************************************************
62 % Constant properties
63 %************************************************************
64 properties(Constant = true, Hidden = true)
65 numAudiometricFreqs = 6;
66 end
67
68 %% **********************************************************
69 % Dependent visable properties - calculated at runtime
70 %************************************************************
71 properties(Dependent = true, Hidden = false)
72 channelBFs %= 250 * 2.^((0:fNmax)'*params.bwOct);
73 numChannels %= numel(channelBFs);
74 aidOPnice %aid output reformatted to be exactly the same dimensions as the input stimulus
75 end
76
77 %% **********************************************************
78 % Dependent invisable properties - calculated at runtime
79 %************************************************************
80 properties(Dependent = true, Hidden = true)
81 TC_dBO_INTERP % Compression threshold in terms of 2nd filter o/p in dB SPL
82 TM_dBO_INTERP % MOC threshold in terms of 2nd filter o/p in dB SPL
83 bwOct_INTERP
84 DRNLb_INTERP %= ( 2e-5 .* 10.^(TCdBO/20)) .^ (1-DRNLc) ;
85 DRNLc_INTERP
86 mainGain_INTERP %Interp'd and in linear units
87
88 ARthresholdPa %= 20e-6*10^(ARthreshold_dB/20);% Pa thresh for triggering AR
89 stimulusINTERNAL %input stimulus in correct format for the Aid algo
90 end
91
92 %% **********************************************************
93 % Protected properties - The user never needs to set
94 %************************************************************
95 properties(Access = protected)
96 aidOP
97 emlc_z
98
99 %--------------------------------------------------------------
100 % ENUMERATIONS USED IN THE FRAME PROCESSOR
101 %--------------------------------------------------------------
102 enumC_ARb = 0;
103 enumC_ARa = 2;
104 enumC_MOCb = 4;
105 enumC_MOCa = 6;
106
107 % enumC_BPb1 = 8;
108 % enumC_BPa1 = 13;
109 % enumC_BPb2 = 18;
110 % enumC_BPa2 = 23;
111 % enumC_BPb3 = 28;
112 % enumC_BPa3 = 33;
113 % enumC_BPb4 = 38;
114 % enumC_BPa4 = 43;
115
116 enumS_AR = 0;
117
118 % enumS_MOC1 = 1;
119 % enumS_BPin_1_1 = 2;
120 % enumS_BPin_2_1 = 6;
121 % enumS_BPout_1_1 = 10;
122 % enumS_BPout_2_1 = 14;
123 %
124 % enumS_MOC2 = 18;
125 % enumS_BPin_1_2 = 19;
126 % enumS_BPin_2_2 = 23;
127 % enumS_BPout_1_2 = 27;
128 % enumS_BPout_2_2 = 31;
129 % ...
130 end
131
132 %% **********************************************************
133 % methods _ _ _
134 % | | | | | |
135 % _ __ ___ ___| |_| |__ ___ __| |___
136 %| '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __|
137 %| | | | | | __/ |_| | | | (_) | (_| \__ \
138 %|_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/
139 %************************************************************
140
141 methods
142 %% **********************************************************
143 % Constructor
144 %************************************************************
145 function obj = EssexAid_WrapClass(sr, stimulus)
146
147 if nargin > 0
148 obj.sr = sr;
149 end
150
151 if nargin > 1
152 obj.stimulusUSER = stimulus;
153 else
154 obj.stimulusUSER = obj.pipSequence(obj.sr);
155 end
156 end
157
158 %% **********************************************************
159 % Get method for channelBFs
160 %************************************************************
161 function value = get.channelBFs(obj)
162 fNmax = 5/obj.bwOct;
163 value = 250 * 2.^((0:fNmax)'*obj.bwOct);
164 end
165
166 %% **********************************************************
167 % Get method for numChannels
168 %************************************************************
169 function value = get.numChannels(obj)
170 value = numel(obj.channelBFs);
171 end
172
173 %% **********************************************************
174 % Get method for ARthresholdPa
175 %************************************************************
176 function value = get.ARthresholdPa(obj)
177 value = 20e-6*10^(obj.ARthreshold_dB/20);% Pa thresh for triggering AR
178 end
179
180 %% **********************************************************
181 % Get method for TC_dBO_INTERP
182 %************************************************************
183 function value = get.TC_dBO_INTERP(obj)
184 TC_dBO = obj.audiometry_dB - obj.mainGain_dB + obj.TC_dBHL;
185 value = obj.interpPars(TC_dBO, obj.numChannels);
186 end
187
188 %% **********************************************************
189 % Get method for TM_dBO_INTERP
190 %************************************************************
191 function value = get.TM_dBO_INTERP(obj)
192 TM_dBO = obj.audiometry_dB - obj.mainGain_dB + obj.TM_dBHL;
193 value = obj.interpPars(TM_dBO, obj.numChannels);
194 end
195
196 %% **********************************************************
197 % Get method for bwOct_INTERP
198 %************************************************************
199 function value = get.bwOct_INTERP(obj)
200 value = repmat(obj.bwOct, 1, obj.numChannels);
201 end
202
203 %% **********************************************************
204 % Get method for DRNLb_INTERP
205 %************************************************************
206 function value = get.DRNLb_INTERP(obj)
207 value = ( 2e-5 .* 10.^(obj.TC_dBO_INTERP/20)) .^ (1-obj.DRNLc_INTERP);
208 end
209
210 %% **********************************************************
211 % Get method for DRNLc_INTERP
212 %************************************************************
213 function value = get.DRNLc_INTERP(obj)
214 value = obj.interpPars(obj.DRNLc, obj.numChannels);
215 end
216
217 %% **********************************************************
218 % Get method for mainGain_INTERP
219 %************************************************************
220 function value = get.mainGain_INTERP(obj)
221 mainGainLin = 10.^(obj.mainGain_dB/20); %lin units
222 value = obj.interpPars(mainGainLin, obj.numChannels);
223 end
224
225 %% ***********************************************************
226 % Get method for stimulus
227 % -----------------------
228 % The hearing aid expects a stereo signal, as the MOC control is
229 % linked for left and right channels. It would be more efficient to
230 % use a mono version of the aid for simulation in Matlab. However,
231 % I always want to use the exact same code for the hardware in the
232 % lab and current simulations. This code will make a mono signal
233 % stereo if needs be and/or rotate to 2xN array.
234 %*************************************************************
235 function value = get.stimulusINTERNAL(obj)
236 [nRows, nCols] = size(obj.stimulusUSER);
237
238 % Assume that the stimulus duration is greater than 2 samples.
239 % Therefore the number of channels is the min dim.
240 [nChans, I] = min([nRows nCols]);
241
242 if nChans == 2
243 if I == 2
244 value = obj.stimulusUSER;
245 else
246 value = obj.stimulusUSER';
247 end
248 elseif nChans == 1 %Just to be explicit
249 if I == 2
250 value = [obj.stimulusUSER obj.stimulusUSER];
251 else
252 value = [obj.stimulusUSER; obj.stimulusUSER]';
253 end
254 end
255 end
256
257 %% ***********************************************************
258 % Get method for aid output
259 % -----------------------
260 % This get method is linked to the above internal stimulus method
261 % and allows the user to extract the hearing aid output in exactly
262 % the same shape and size as the original input stimulus. This is
263 % very useful for the speech recognition work and presumably
264 % for multithreshold also.
265 %*************************************************************
266 function value = get.aidOPnice(obj)
267 if ~isempty(obj.aidOP)
268 [nRows, nCols] = size(obj.stimulusUSER);
269
270 % Assume that the stimulus duration is greater than 2 samples.
271 % Therefore the number of channels is the min dim.
272 [nChans, I] = min([nRows nCols]);
273
274 %** The aid output will ALWAYS be a 2xN array **
275 %The fist job is to remove trailing zeros that may have been
276 %introduced by the framing process
277 aidOPtruncated = obj.aidOP(:, 1:max([nRows nCols]));
278
279 %The next task is to arrange the op like the ip
280 if nChans == 2
281 if I == 1
282 value = aidOPtruncated;
283 else
284 value = aidOPtruncated';
285 end
286 elseif nChans == 1 %Just to be explicit
287 if I == 1
288 value = aidOPtruncated(1,:);
289 else
290 value = aidOPtruncated(1,:)';
291 end
292 end
293 else % ---- of if isempty statement
294 value = [];
295 end
296 end
297
298 %% ***********************************************************
299 % *** Set methods ***
300 % -----------------------
301 % This is a bunch of unexciting error hunting functions. They also
302 % flush the aid output if any parameters change. Therefore,
303 % processStim will have to be called explicity by the user once
304 % again.
305 %*************************************************************
306 function obj = set.stimulusUSER(obj,value)
307 [nRows, nCols] = size(value);
308
309 % Assume that the stimulus duration is greater than 2 samples.
310 % Therefore the number of channels is the min dim.
311 nChans = min([nRows nCols]);
312 assert(nChans<3 && nChans, 'Number of stimulus channels must be 1 or 2')
313
314 obj = obj.flushAidData; %flush any previous hearing aid data if the input stimulus changes
315 obj.stimulusUSER = value;
316 end
317 function obj = set.sr(obj,value)
318 assert(value>=20e3 && value<=192e3, 'sr must be between 20 and 192 kHz')
319 obj = obj.flushAidData;
320 obj.sr = value;
321 end
322 function obj = set.numSamples(obj,value)
323 assert(value>=48 && value<=6912, 'must be between 48 and 6912 samples')
324 obj = obj.flushAidData;
325 obj.numSamples = value;
326 end
327 function obj = set.audiometry_dB(obj,value)
328 [nRows,nCols] = size(value);
329 assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>
330 obj = obj.flushAidData;
331 obj.audiometry_dB = value;
332 end
333 function obj = set.mainGain_dB(obj,value)
334 [nRows,nCols] = size(value);
335 assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>
336 obj = obj.flushAidData;
337 obj.mainGain_dB = value;
338 end
339 function obj = set.TC_dBHL(obj,value)
340 [nRows,nCols] = size(value);
341 assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>
342 obj = obj.flushAidData;
343 obj.TC_dBHL = value;
344 end
345 function obj = set.TM_dBHL(obj,value)
346 [nRows,nCols] = size(value);
347 assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>
348 obj = obj.flushAidData;
349 obj.TM_dBHL = value;
350 end
351 function obj = set.DRNLc(obj,value)
352 [nRows,nCols] = size(value);
353 assert(nRows==obj.numAudiometricFreqs && nCols==1, 'must be 6x1 column vector') %#ok<MCSUP>
354 assert(all(value)>=0 && all(value)<=1, 'all DRNLc values must be between 0 and 1')
355 obj = obj.flushAidData;
356 obj.DRNLc = value;
357 end
358 function obj = set.ARtau(obj,value)
359 assert(value>=1e-3 && value<=1, 'must be between 1e-3 and 1s')
360 obj = obj.flushAidData;
361 obj.ARtau = value;
362 end
363 function obj = set.ARthreshold_dB(obj,value)
364 assert(value>0, 'set AR to a high value to disable it')
365 obj = obj.flushAidData;
366 obj.ARthreshold_dB = value;
367 end
368 function obj = set.MOCtau(obj,value)
369 assert(value>=1e-3 && value<=2, 'must be between 1e-3 and 2s')
370 obj = obj.flushAidData;
371 obj.MOCtau = value;
372 end
373 function obj = set.MOCfactor(obj,value)
374 assert(value>=0 && value<=1, 'must be between 0 and 1')
375 obj = obj.flushAidData;
376 obj.MOCfactor = value;
377 end
378 function obj = set.bwOct(obj,value)
379 assert(value==1/1 || value==1/2 || value==1/3 || value==1/4 || value==1/5, 'must be one of 1./(1:5)')
380 obj = obj.flushAidData;
381 obj.bwOct = value;
382 end
383 function obj = set.filterOrder(obj,value)
384 assert(value>0 && value<5, 'must be one of 1:4')
385 obj = obj.flushAidData;
386 obj.filterOrder = value;
387 end
388 function obj = set.useGTF(obj,value)
389 obj = obj.flushAidData;
390 obj.useGTF = value;
391 end
392
393 %% **********************************************************
394 % flushAidData
395 % This second function is a workaround allowing a set method to
396 % change another property value.
397 %************************************************************
398 function obj = flushAidData(obj)
399 obj.aidOP = [];
400 obj.MOCrecord = [];
401 end
402
403
404 %% **********************************************************
405 % OVERLOADED plot method
406 %************************************************************
407 function plot(obj)
408 clf
409 sig2dBSPL = @(sig)20*log10(abs(sig/20e-6)+(1/(2^32)));
410 dt = 1/obj.sr;
411 tAxis = dt:dt:dt*size(obj.stimulusINTERNAL,1);
412
413 subplot(2,1,1)
414 plot(tAxis(1:length(obj.stimulusUSER)), sig2dBSPL(obj.stimulusUSER), 'k')
415 if ~isempty(obj.aidOPnice)
416 hold on
417 plot(tAxis(1:length(obj.stimulusUSER)), sig2dBSPL(obj.aidOPnice), 'r')
418 end
419 ylim([0 100])
420 xlim([0 tAxis(length(obj.stimulusUSER))])
421 title('Level response')
422 xlabel('Time in seconds')
423 ylabel('Level in dB SPL')
424
425 subplot(2,1,2)
426 if ~isempty(obj.MOCrecord)
427 imagesc(tAxis, 1:obj.numChannels, flipud(-20*log10(obj.MOCrecord)))
428 colorbar
429 end
430 title('MOC attenuation')
431 xlabel('Time in seconds')
432 ylabel('Band frequency in Hz')
433 numSpacers = 1 + (obj.numChannels-numel(obj.DRNLc)) / (numel(obj.DRNLc)-1);
434 set(gca, 'YTick', 1:numSpacers:obj.numChannels);
435 set(gca, 'YTickLabel', num2str(flipud([250; 500; 1000; 2000; 4000; 8000])));
436 end% ------ OVERLOADED plot method
437
438 %% **********************************************************
439 % OVERLOADED soundsc method
440 %************************************************************
441 function soundsc(obj)
442 soundsc(obj.aidOPnice, obj.sr)
443 end
444
445 %% **********************************************************
446 % processStim
447 %************************************************************
448 function obj = processStim(obj)
449 %--------------------------------------------------------------
450 % EMULATION OF THE GUI PARAMETER CONVERSIONS
451 %--------------------------------------------------------------
452 biggestNumSamples = obj.numSamples;
453
454 filterStatesL = (zeros(3000,1));
455 filterStatesR = filterStatesL;
456 filterCoeffs = (zeros(5000,1));
457
458 %filter coefficients
459 ARcutOff=1/(2*pi*obj.ARtau);
460 [b,a] = butter(1,ARcutOff/(obj.sr/2));
461 filterCoeffs(obj.enumC_ARb+1:obj.enumC_ARb+2) = b;
462 filterCoeffs(obj.enumC_ARa+1:obj.enumC_ARa+2) = a;
463
464 MOCcutOff=1/(2*pi*obj.MOCtau);
465 [bMOC,aMOC] = butter(1,MOCcutOff/(obj.sr/2));
466 filterCoeffs(obj.enumC_MOCb+1:obj.enumC_MOCb+2) = bMOC;
467 filterCoeffs(obj.enumC_MOCa+1:obj.enumC_MOCa+2) = aMOC;
468
469
470 for filterCount = 1:obj.numChannels
471 %-----------------------------------
472 % nonlinear path - filter bws
473 %-----------------------------------
474 lowerCutOff=obj.channelBFs(filterCount)*2^(-obj.bwOct_INTERP(filterCount)/2);
475 upperCutOff=obj.channelBFs(filterCount)*2^( obj.bwOct_INTERP(filterCount)/2);
476
477 if obj.useGTF
478 bwHz = upperCutOff - lowerCutOff;
479 [b_DRNL,a_DRNL] = obj.gammatone(bwHz, obj.channelBFs(filterCount), 1/obj.sr);
480 filterCoeffs(10*(filterCount-1)+9 :10*(filterCount-1)+10) = b_DRNL;
481 filterCoeffs(10*(filterCount-1)+14:10*(filterCount-1)+16) = a_DRNL;
482 else
483 [b_DRNL,a_DRNL] = butter(2,[lowerCutOff upperCutOff]/(obj.sr/2));
484 filterCoeffs(10*(filterCount-1)+9 :10*(filterCount-1)+13) = b_DRNL;
485 filterCoeffs(10*(filterCount-1)+14:10*(filterCount-1)+18) = a_DRNL;
486 end
487 end
488
489 %--------------------------------------------------------------
490 % EMULATION OF THE IO CALLBACK THREAD
491 %--------------------------------------------------------------
492 frameBufferL = buffer(obj.stimulusINTERNAL(:,1), obj.numSamples);
493 frameBufferR = buffer(obj.stimulusINTERNAL(:,2), obj.numSamples);
494 nFrames = size(frameBufferL,2);
495
496 pad = zeros(1,biggestNumSamples-obj.numSamples);
497 ARampL=ones(1,biggestNumSamples);
498 ARampR = ARampL;
499 MOCcontrol = ones(obj.numChannels, biggestNumSamples);
500
501 peakIPL = zeros(5,1);
502 peakOPL = peakIPL;
503 rmsIPL = peakIPL;
504 rmsOPL = peakIPL;
505
506 peakIPR = peakIPL;
507 peakOPR = peakIPL;
508 rmsIPR = peakIPL;
509 rmsOPR = peakIPL;
510
511 MOCend = zeros(obj.numChannels,1);
512
513 op = [];
514 moc= [];
515 for nn = 1:nFrames
516 frameBufferPadL = [frameBufferL(:,nn)' pad];
517 frameBufferPadR = [frameBufferR(:,nn)' pad];
518
519 [ outBufferL, outBufferR, filterStatesL, filterStatesR, ARampL, ARampR, MOCend, peakIPL, peakOPL, rmsIPL, rmsOPL, peakIPR, peakOPR, rmsIPR, rmsOPR, MOCcontrol ] =...
520 EssexAidProcessVFrameSwitchable( ...
521 frameBufferPadL,...
522 frameBufferPadR,...
523 filterStatesL,...
524 filterStatesR,...
525 filterCoeffs,...
526 obj.numChannels,...
527 obj.numSamples,...
528 ARampL,...
529 ARampR,...
530 obj.ARthresholdPa,...
531 obj.filterOrder,...
532 obj.DRNLb_INTERP,...
533 obj.DRNLc_INTERP,...
534 obj.TM_dBO_INTERP,...
535 obj.MOCfactor,...
536 peakIPL,...
537 peakOPL,...
538 rmsIPL,...
539 rmsOPL,...
540 peakIPR,...
541 peakOPR,...
542 rmsIPR,...
543 rmsOPR,...
544 MOCend,...
545 MOCcontrol,...
546 obj.mainGain_INTERP,...
547 obj.useGTF);
548
549
550 outBuffer = ( [outBufferL(:, 1:obj.numSamples); outBufferR(:, 1:obj.numSamples)] );
551 op = [op outBuffer]; %#ok<AGROW>
552 moc= [moc MOCcontrol]; %#ok<AGROW>
553
554 end %End of frame processing emulation loop
555 obj.aidOP = op;
556 obj.MOCrecord=moc;
557
558
559 end %End of process stim method
560
561 end %End of methods block
562
563 %% *********************************************************
564 % _ _ _ _ _ _
565 % | | | | (_) | | | | | |
566 % ___| |_ __ _| |_ _ ___ _ __ ___ ___| |_| |__ ___ __| |___
567 % / __| __/ _` | __| |/ __| | '_ ` _ \ / _ \ __| '_ \ / _ \ / _` / __|
568 % \__ \ || (_| | |_| | (__ | | | | | | __/ |_| | | | (_) | (_| \__ \
569 % |___/\__\__,_|\__|_|\___| |_| |_| |_|\___|\__|_| |_|\___/ \__,_|___/
570 %************************************************************
571
572 methods(Static)
573 %% ********************************************************
574 % pipOut - sequence of tone pips at various levels
575 %**********************************************************
576 function pipOut = pipSequence(sampleRate, freq, dBlevs, pulseDur, silDur)
577 if nargin < 5
578 silDur = 0.3;
579 end
580 if nargin < 4
581 pulseDur = 0.1;
582 end
583 if nargin < 3
584 dBlevs = 20:20:100;
585 end
586 if nargin < 2
587 freq = 500;
588 end
589 if nargin < 1
590 sampleRate = 48e3;
591 end
592
593 dt = 1/sampleRate;
594 tAxis = dt:dt:pulseDur;
595 sPulse = sin(2*pi*freq*tAxis);
596 sPulse = sPulse./sqrt(mean(sPulse.^2));
597 rms2dBspl = @(dBspl)20e-6*10^(dBspl/20); %sneaky short-hand function by (ab)using function handles
598 zPad = zeros(1,ceil(sampleRate*silDur));
599
600 pipOut = [];
601 for nn = 1:numel(dBlevs)
602 pipOut = [ pipOut sPulse*rms2dBspl(dBlevs(nn)) zPad]; %#ok<AGROW>
603 end
604
605 end% ------ OF pipSequence
606
607 %% ********************************************************
608 % interpPars - Linear interpolation of given parameter to mimic GUI
609 % fitting functionality.
610 %**********************************************************
611 function fullArray = interpPars(shortArray, numBands)
612 nGUIbands = numel(shortArray);
613 if numBands == nGUIbands
614 fullArray = shortArray;
615 else
616 numSpacers = (numBands-nGUIbands) / (nGUIbands-1);
617 fullArray = shortArray(1);
618 for nn = 2:nGUIbands
619 fullArray = [fullArray,...
620 repmat(mean([shortArray(nn) shortArray(nn-1)]),1,numSpacers),...
621 shortArray(nn)]; %#ok<AGROW>
622 end
623 end
624 end% ----- OF interpPars
625
626 %% ********************************************************
627 % gammatone - get filter coefficients
628 %**********************************************************
629 function [b,a] = gammatone(bw, cf, dt)
630 phi = 2 * pi * bw * dt;
631 theta = 2 * pi * cf * dt;
632 cos_theta = cos(theta);
633 sin_theta = sin(theta);
634 alpha = -exp(-phi) * cos_theta;
635 b0 = 1.0;
636 b1 = 2 * alpha;
637 b2 = exp(-2 * phi);
638 z1 = (1 + alpha * cos_theta) - (alpha * sin_theta) * 1i;
639 z2 = (1 + b1 * cos_theta) - (b1 * sin_theta) * 1i;
640 z3 = (b2 * cos(2 * theta)) - (b2 * sin(2 * theta)) * 1i;
641 tf = (z2 + z3) / z1;
642 a0 = abs(tf);
643 a1 = alpha * a0;
644
645 a = [b0, b1, b2];
646 b = [a0, a1];
647 end% ------ OF gammatone
648 end% ------ OF static methods
649
650 end %End of classdef
651