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