diff aim-mat/modules/pcp/gm2002/ff_design.m @ 0:74dedb26614d

Initial checkin of AIM-MAT version 1.5 (6.4.2011).

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/aim-mat/modules/pcp/gm2002/ff_design.m	Fri May 20 12:32:31 2011 +0100
@@ -0,0 +1,176 @@
+function fir_eq=genstonefir(fs,type)
+% adapted version von S.Bleeck
+%
+%%% May 2001 M.A.Stone.  To design FIR filters for fir_corr.c
+%%% produces output on screen, and to file dummy.??k where ?? is clock freq in kHz
+%%% select variables below, such as ntaps (output is ntaps +1) fs,
+%%% and inverse =1 sets for inverse filter, 0 for normal filter.
+%%% and eq characteristic is set by choosing appropriate variables:
+%%%% function_str and dB corrn.  Ff_ed , Df_ed, Midear, and Hz come
+%%%% from a separate file [all_corrns.m]. Insert new characteristics
+%%%% in there, eg ITU_erp_drp/ITU_Hz
+
+%%%  FORMAT for output file
+%%%  first line is comment
+%%%% second line is number of taps (preferably odd)
+%%%% third and subsequent lines are filter taps, one per line, floating point
+
+all_corrns; %%%%%% external file for reference corrections, hz, midear, ff_ed, diffuse
+%%%%%% NB, midear response has limit/flatten-off at lowest freqs to prevent enormous changes < 25 Hz
+%%%%%%%% design parameters here
+%%% NB sometomes for inverse, cannot have ntaps too high: claims index error in fir2.
+% fs = 50000;  %%%% sampling rate
+ntaps = 1+2*(round(fs/24)); %% always odd
+nFFT = 2.^(nextpow2(ntaps) + 1); %% FIR size is ntaps + 1, otherwise delay has extra half-sample
+%% more taps require kaiser beta to be higher
+inverse = 0; %% options 0/1 : whether to do forward or inverse filter
+
+posh_print = 0; %%% just if we want publication figure, so no output file
+
+beta = 6; %%% used to window INVERSE filter shape, and reduce ripple
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%% Uncomment which of the three sections below you require
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+switch type
+	case 1
+		%%%%%%% Section 1
+		function_str= sprintf(' Frontal free-field to cochlea correction,fs=%d',fs);
+		if fs <= 32000
+			dBcorrn = Ff_ed - Midear; %%%%% result in dB 
+		else %% for high fs use truer version, without low freq fiddle
+			dBcorrn = Ff_ed - MidearAES; %%%%% result in dB 
+		end
+	case 2
+		%%%%%%% Section 2
+		function_str= sprintf(' Diffuse-field to cochlea correction,fs=%d',fs);
+		if fs <= 32000
+			dBcorrn = Diffuse - Midear; %%%%% result in dB  NB midear is inverted !!
+		else %% for high fs use truer version, without low freq fiddle
+			dBcorrn = Diffuse - MidearAES; %%%%% result in dB
+		end
+		
+	case 3
+		%%%%%% Section 3
+		%%%%% ITU corrections for telephony.
+		function_str= sprintf(' ITU Ear Ref Pnt via Drum Ref Pnt to cochlea,fs=%d',fs);
+		ITU_on_Hz = interp1(ITU_Hz,ITU_erp_drp,Hz,'spline');  %%%% corrn on linear frequency spacing
+		dBcorrn = ITU_on_Hz - Midear; %%%%% result in dB.  NB midear is inverted !!
+end
+%%%% END OF VARIABLE USER ENTRY/SET-UP
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%% rest is just calculations
+if ~inverse
+    filename = sprintf('ff.%dk',floor(fs/1000));
+else
+    filename = sprintf('iff.%dk',floor(fs/1000));
+end
+% if ~posh_print
+% 	txt_file = fopen(filename,'wt');
+% end
+%%%% spacing for linear frequency design grid
+deltaf  = (fs)/nFFT;
+
+linf = 0:deltaf:fs/2; %%%% frequencies of FFT bins (DC-nyq-1)
+
+if inverse == 1
+	dBcorrn = -dBcorrn;
+	title_str = strcat('INVERSE',function_str);   
+else
+   title_str= function_str;
+   title_str=' ';
+end
+
+if Hz(end) < fs/2  %% handle interpolation of high fs data
+    dBcorrn_linf = interp1([Hz fs/2],[dBcorrn dBcorrn(end)],linf,'linear');  %%%% corrn on linear frequency spacing
+else
+    dBcorrn_linf = interp1(Hz,dBcorrn,linf,'linear');  %%%% corrn on linear frequency spacing
+end
+%%%%%% limit/flatten-off at lowest freqs
+dBcorrn_linf_orig = dBcorrn_linf;
+%%%dBcorrn_linf(2) = dBcorrn_linf(3) - (dBcorrn_linf(3)-dBcorrn_linf(2))/2;
+%%%dBcorrn_linf(1) = dBcorrn_linf(2);
+
+[smth_b smth_a] = butter(4,.5); %% smooth to control roughness
+eq_linf = filtfilt(smth_b,smth_a,dBcorrn_linf);
+
+% if posh_print, figure(1); plot(linf,eq_linf,'r','linewidth',1.8); hold on; end
+
+%%%%% design fir filter: NB taming of response by (gentle) Kaiser window
+if inverse
+    halfwid = 10.^(eq_linf/20.);
+    npi = pi*mod((0:nFFT/2),2);  %% include phase shift to put response in middle of aperture
+    phase_shift = exp(i*npi);
+    halfwid = halfwid .* phase_shift;
+    fullwid = [halfwid conj(halfwid(nFFT/2:-1:2))];
+    t_filt = real(ifft(fullwid));
+    %%%%%%%%figure(2); plot(real(t_filt)); figure(1);
+    ntaps2 = floor(ntaps/2);  %% extract relevant portion
+    fir_eq = t_filt((nFFT/2+1)-ntaps2 : (nFFT/2+1)+ntaps2);
+    fir_eq = fir_eq.*kaiser(ntaps,beta)';
+else
+    %TCW AIM 2006 naps -> ntaps+1 to remove error from fir2 which makes
+    %the filter order even anyway. (Ntaps is always odd, see line 23).
+    fir_eq = fir2(ntaps+1,linf./(fs/2),10.^(eq_linf/20.)); % f= 1 is Nyquist
+end
+    
+
+%%%%% and plot response
+if posh_print ~=0
+	[hz,fz] = freqz(fir_eq,1,16384,fs);
+
+%  GET AN OUTPUT FILE
+%    outfile=fopen('frq_res.ff','w');
+%    fprintf(outfile,'%.4f,%.4f\n',[fz,20*log10(abs(hz))]');
+%    fclose(outfile);
+	plot(fz,20*log10(abs(hz)),'b','linewidth',1.8);
+	set(gca,'box','on'); %%%% default with R12 is off
+	title(title_str,'fontsize',13); xlabel('Frequency (Hz)','fontsize',11); ylabel('Relative transmission (dB)','fontsize',11);
+	set(gca,'TickDirMode','manual','TickLength',[0 0]); %% turn off ticking
+
+	xfl = 20-.1; xfh = fs/2;
+	xticking = [20 50 100 200 500 1000 2000 5000 10000];
+    if fs/2 > 20e3, xticking = [xticking 20e3]; end
+    if fs/2 > 50e3, xticking = [xticking 50e3]; end % no point in any higher
+        
+	set(gca,'xlim',[xfl xfh],'xscale','log');
+	set(gca,'xtickmode','manual','xtick',xticking,'xticklabel',xticking); 
+
+	if inverse, dBl = -10; dBh = 40; else dBl = -40; dBh = 10; end
+	yticking = [dBl:5:dBh];
+	set(gca,'linewidth',1.3,'ylim',[dBl dBh],'fontsize',11);
+	set(gca,'ytickmode','manual','ytick',yticking,'yticklabel',yticking); 
+
+%%	grid on ; set(gca,'GridLineStyle','-');
+%%%%%% to overcome bugs in MATLAB with xscale producing extra ticks 20-06-2001
+	for ix = 1:length(xticking) %% draw ylines
+		line([xticking(ix) xticking(ix)],[min(yticking), max(yticking)],'linewidth',0.6,'linestyle','--');
+	end
+	for ix = 1:length(yticking) %% draw xlines
+		line([xfl xfh],[yticking(ix), yticking(ix)],'linewidth',0.6,'linestyle','--');
+	end
+
+else
+% 	freqz(fir_eq,1,8192,fs);
+% 	subplot(2,1,1); set(gca,'xlim',[10 fs/2],'xscale','log');
+% 	hold on ; grid on;
+% 	semilogx(linf,dBcorrn_linf_orig,'r');
+% 	title(title_str); xlabel('frequency (Hz)'); ylabel('dB (red=target, blue=actual)');
+% 	subplot(2,1,1); hold off
+% 	subplot(2,1,2); hold off
+end
+
+% hold off  % for figure(1)
+
+% %%%% print out design values to file (and was screen)
+% if posh_print ~=0
+% 	fprintf(1,'\nThis version has also created the file %s\n',filename);
+% 	fprintf(1,'%s\n',function_str);
+% 	fprintf(1,'%d\n',length(fir_eq));
+% %%	fprintf(1,'%11.8f\n',fir_eq);
+% 	fprintf(txt_file,'%s\n',function_str);
+% 	fprintf(txt_file,'%d\n',length(fir_eq));
+% 	fprintf(txt_file,'%f\n',fir_eq);
+% 	fclose(txt_file);
+% end
+% 
+% WriteDSAMFIRParFile(fir_eq, fs, inverse);