tomwalters@0
|
1 %
|
tomwalters@0
|
2 % Produce compensation filter to simulate outer/middle ear
|
tomwalters@0
|
3 % IRINO Toshio
|
tomwalters@0
|
4 % 29 Aug. 1996 (check on 14 May 1997 )
|
tomwalters@0
|
5 % 8 Jan. 2002 (Multiply Win for avoid sprious)
|
tomwalters@0
|
6 % 19 Nov. 2002 (remez : even integer)
|
tomwalters@0
|
7 %
|
tomwalters@0
|
8 % It is a linear phase filter for the ELC/MAF/MAP correction.
|
tomwalters@0
|
9 % see OutMidCrct.m
|
tomwalters@0
|
10 %
|
tomwalters@0
|
11 % function [FIRCoef] = OutMidCrctFilt(StrCrct,SR);
|
tomwalters@0
|
12 % INPUT StrCrct: String for Correction ELC/MAF/MAP
|
tomwalters@0
|
13 % SR: Sampling Rate
|
tomwalters@0
|
14 % SwPlot: SwPlot
|
tomwalters@0
|
15 % OUTPUT FIRCoef: FIR filter coefficients
|
tomwalters@0
|
16 %
|
tomwalters@0
|
17 function [FIRCoef] = OutMidCrctFilt(StrCrct,SR,SwPlot);
|
tomwalters@0
|
18
|
tomwalters@0
|
19 if nargin < 2, help OutMidCrctFilt; end;
|
tomwalters@0
|
20 if nargin < 3, SwPlot = 1; end;
|
tomwalters@0
|
21
|
tomwalters@0
|
22 if length(StrCrct)~=3,
|
tomwalters@0
|
23 error('Specifiy correction in 3 characters: ELC / MAF / MAP.');
|
tomwalters@0
|
24 end;
|
tomwalters@0
|
25 if ~(strcmp(upper(StrCrct(1:3)), 'ELC') | ...
|
tomwalters@0
|
26 strcmp(upper(StrCrct(1:3)),'MAF') ...
|
tomwalters@0
|
27 | strcmp(upper(StrCrct(1:3)),'MAP')),
|
tomwalters@0
|
28 error('Specifiy correction: ELC / MAF / MAP.');
|
tomwalters@0
|
29 end;
|
tomwalters@0
|
30
|
tomwalters@0
|
31 Nint = 1024;
|
tomwalters@0
|
32 % Nint = 0; % No spline interpolation: NG no convergence at remez
|
tomwalters@0
|
33 [crctPwr freq] = OutMidCrct(StrCrct,Nint,SR,0);
|
tomwalters@0
|
34 crct = sqrt(crctPwr);
|
tomwalters@0
|
35
|
tomwalters@0
|
36 %% FIRCoef = remez(50/16000*SR,freq/SR*2,crct); % NG
|
tomwalters@0
|
37 %% FIRCoef = remez(300/16000*SR,freq/SR*2,crct); % Original
|
tomwalters@0
|
38 % FIRCoef = remez(LenCoef/16000*SR,freq/SR*2,crct); % when odd num : warning
|
tomwalters@0
|
39 %% modified on 8 Jan 2002, 19 Nov 2002
|
tomwalters@0
|
40 LenCoef = 200; % ( -45 dB) <- 300 (-55 dB)
|
tomwalters@0
|
41 FIRCoef = remez(fix(LenCoef/16000*SR/2)*2,freq/SR*2,crct); % even number only
|
tomwalters@0
|
42 Win = TaperWindow(length(FIRCoef),'han',LenCoef/10);
|
tomwalters@0
|
43 % Necessary to avoid sprious
|
tomwalters@0
|
44 FIRCoef = Win.*FIRCoef;
|
tomwalters@0
|
45
|
tomwalters@0
|
46 if SwPlot==1
|
tomwalters@0
|
47 [frsp freq2] = freqz(FIRCoef,1,Nint,SR);
|
tomwalters@0
|
48 subplot(2,1,1)
|
tomwalters@0
|
49 plot(FIRCoef);
|
tomwalters@0
|
50 subplot(2,1,2)
|
tomwalters@0
|
51 plot(freq2,abs(frsp),freq,crct,'--')
|
tomwalters@0
|
52 % plot(freq2,20*log10(abs(frsp)),freq,20*log10(crct))
|
tomwalters@0
|
53
|
tomwalters@0
|
54 ELCError = mean((abs(frsp) - crct).^2)/mean(crct.^2);
|
tomwalters@0
|
55 ELCErrordB = 10*log10(ELCError) % corrected
|
tomwalters@0
|
56 if ELCErrordB > -30,
|
tomwalters@0
|
57 disp(['Warning: Error in ELC correction = ' ...
|
tomwalters@0
|
58 num2str(ELCErrordB) ' dB > -30 dB'])
|
tomwalters@0
|
59 end;
|
tomwalters@0
|
60 end;
|
tomwalters@0
|
61
|
tomwalters@0
|
62 return;
|