function per = aud_test()
%
% auditory perception test program
%  gomi@cslab 2000.8.2  
%

f_t = 250*2.^((0:1:11)/3)
a_t = 0.005:0.01:0.3
a_t = exp(-5.4:0.3:-1.2)

per = zeros(size(f_t,2),size(a_t,2));
for i = 1:size(f_t,2)
 for j = 1:size(a_t,2)
   per(i,j) = 2;
   while per(i,j) == 2
     aud_test_sub(f_t(i),a_t(j));
     per(i,j) = input ('continuous tone ? (y:1, n:0, repeat:2) -> ');
   end
   if per(i,j) == 0
     break
   end
 end
end

[o_f, n_f] = aud_test_sub(f_t(1),0);

dB = 20*log10((a_t/sqrt(2))/sqrt(mean(n_f.^2)));
dBmat = (ones(size(per,1),1)*dB).*per;
for i =1:size(per,1)
  if ~isempty(max(find(per(i,:)==1))) 
    dBmax(i) = dBmat(i, max(find(per(i,:)==1)));
  end
end

plot(f_t(find(dBmax<0)),dBmax(find(dBmax<0)));
set(gca,'xlim',[200,4000])
set(gca,'xtick',[200,400,800,1000,2000,4000])
xlabel('freq [Hz]')
ylabel('RMS tone/noise [dB]')


%%%%%%%%%%%%%%%%%%%%%%
function [o_f,n_f,s] = aud_test_sub(f,a)

% a: tone amplitude
% f: tone frequency

f_s = 8000;  % Hz sampling freq

if 0
a = 0.3     % tone amplitude
f = 1000    % Hz  tone frequency
end
fn = 1000;    % Hz noise freq center
f
len = 2.0;   % sec   sound length
n_dur = 0.2; % sec   noise duration
s_dur = 0.2; % sec   tone duration

f_sH = f_s/2;

fil_sl = f*2^(-1/24);  % tone filter cutoff freq low
fil_sh = f*2^(1/24);   % tone filter cutoff freq high

fil_nl = fn*2^(-1/6);  % noise filter cutoff freq low
fil_nh = fn*2^(1/6);   % noise filter cutoff freq high

Forder = 4; % filter order 

%[fil_sB,fil_sA] = butter(Forder,[fil_sl,fil_sh]/f_sH);
[fil_nB,fil_nA] = butter(Forder,[fil_nl,fil_nh]/f_sH);


mask = ones(1,len*f_s);
mask_c =  cos(pi*100*[0:1/f_s:(0.01-1/f_s)])/2+0.5;
mask_cp = 0.01*f_s;


for i = 1:floor(len/(n_dur+s_dur))
  if ((s_dur+n_dur)*i*f_s-1+mask_cp) < size(mask,2)
    mask(round(((s_dur)+(s_dur+n_dur)*(i-1))*f_s): ...
	 round((s_dur+n_dur)*i*f_s-1+mask_cp))=0;
    mask(round(((s_dur)+(s_dur+n_dur)*(i-1))*f_s): ...
	 round(((s_dur)+(s_dur+n_dur)*(i-1))*f_s)+mask_cp-1)=mask_c;
    mask(round((s_dur+n_dur)*i*f_s): ...
	 round((s_dur+n_dur)*i*f_s-1+mask_cp))=-mask_c+1;
  else
    mask(round(((s_dur)+(s_dur+n_dur)*(i-1))*f_s): ...
	 round((s_dur+n_dur)*i*f_s-1))=0;
    mask(round(((s_dur)+(s_dur+n_dur)*(i-1))*f_s): ...
	 round(((s_dur)+(s_dur+n_dur)*(i-1))*f_s)+mask_cp-1)=mask_c;
  end
end


s =  a*sin(2*pi*f*[0:1/f_s:(len-1/f_s)]);
sm= (s.*mask);
%sm_f =  filter(fil_sB,fil_sA,sm);
sm_f =  sm;

%sound(sm_f,f_s)

n = randn(size(s));
n_f =  filter(fil_nB,fil_nA,n);
nm_f = n_f.*(mask*(-1)+1);


%sound(nm_f,f_s)

o_f=sm_f+nm_f;
sound(o_f,f_s)

if 0
ffn=fft(nm_f);
ffs=fft(s);
logplot(abs(ffn(1,1:4000)))
semilogy((1:8000)/2,abs(ffn(1,1:8000)))
semilogy(abs(ffn))
hold on
semilogy((1:8000)/2,abs(ffs(1,1:8000)),'r')
end
%%%%%%%%%%%%%%%%%%%%%%