tomwalters@0: %
tomwalters@0: %	Produce compensation filter to simulate outer/middle ear
tomwalters@0: %	IRINO Toshio
tomwalters@0: %	29 Aug. 1996 (check on 14 May  1997 )
tomwalters@0: %	8 Jan. 2002  (Multiply Win for avoid sprious)
tomwalters@0: %	19 Nov. 2002 (remez : even integer)
tomwalters@0: %
tomwalters@0: %	It is a linear phase filter for the ELC/MAF/MAP correction.
tomwalters@0: %	see OutMidCrct.m 
tomwalters@0: %
tomwalters@0: %	function [FIRCoef] = OutMidCrctFilt(StrCrct,SR);
tomwalters@0: %	INPUT	StrCrct: String for Correction ELC/MAF/MAP
tomwalters@0: %		SR: 	 Sampling Rate
tomwalters@0: %		SwPlot:  SwPlot
tomwalters@0: %	OUTPUT  FIRCoef: FIR filter coefficients 
tomwalters@0: %
tomwalters@0: function [FIRCoef] = OutMidCrctFilt(StrCrct,SR,SwPlot);
tomwalters@0: 
tomwalters@0: if nargin < 2, help OutMidCrctFilt; end;
tomwalters@0: if nargin < 3, SwPlot = 1; end;
tomwalters@0: 
tomwalters@0: if length(StrCrct)~=3, 
tomwalters@0: 	error('Specifiy correction in 3 characters: ELC / MAF / MAP.'); 
tomwalters@0: end;
tomwalters@0: if ~(strcmp(upper(StrCrct(1:3)), 'ELC') | ...
tomwalters@0: 	 strcmp(upper(StrCrct(1:3)),'MAF') ...
tomwalters@0: 	| strcmp(upper(StrCrct(1:3)),'MAP')),
tomwalters@0: 	error('Specifiy correction: ELC / MAF / MAP.'); 
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: Nint = 1024;
tomwalters@0: % Nint = 0; % No spline interpolation:  NG no convergence at remez
tomwalters@0: [crctPwr freq] = OutMidCrct(StrCrct,Nint,SR,0);
tomwalters@0: crct = sqrt(crctPwr);
tomwalters@0: 
tomwalters@0: %% FIRCoef = remez(50/16000*SR,freq/SR*2,crct); % NG 
tomwalters@0: %% FIRCoef = remez(300/16000*SR,freq/SR*2,crct); % Original
tomwalters@0: % FIRCoef = remez(LenCoef/16000*SR,freq/SR*2,crct); % when odd num : warning
tomwalters@0: %% modified on 8 Jan 2002, 19 Nov 2002
tomwalters@0: LenCoef = 200; %  ( -45 dB) <- 300 (-55 dB)
tomwalters@0: FIRCoef = remez(fix(LenCoef/16000*SR/2)*2,freq/SR*2,crct);  % even number only
tomwalters@0: Win     = TaperWindow(length(FIRCoef),'han',LenCoef/10); 
tomwalters@0: 	    % Necessary to avoid sprious
tomwalters@0: FIRCoef = Win.*FIRCoef;
tomwalters@0: 
tomwalters@0: if SwPlot==1
tomwalters@0: 	[frsp freq2] = freqz(FIRCoef,1,Nint,SR);
tomwalters@0: 	subplot(2,1,1)
tomwalters@0: 	plot(FIRCoef);
tomwalters@0: 	subplot(2,1,2)
tomwalters@0: 	plot(freq2,abs(frsp),freq,crct,'--')
tomwalters@0: 	%	plot(freq2,20*log10(abs(frsp)),freq,20*log10(crct))
tomwalters@0: 	
tomwalters@0: 	ELCError = mean((abs(frsp) - crct).^2)/mean(crct.^2);
tomwalters@0: 	ELCErrordB = 10*log10(ELCError)          % corrected 
tomwalters@0: 	if ELCErrordB > -30,
tomwalters@0: 	    disp(['Warning: Error in ELC correction = ' ...
tomwalters@0: 		    num2str(ELCErrordB) ' dB > -30 dB'])
tomwalters@0: 	end;
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: return;