% This is the code we used to generate notched noise that will elicit
% Zwicker tone in some subjects. Our data shows that about 1/3 of normal
% subjects will hear the Zwicker tone for at least one of the three
% band-gaps presented here. Subjects that report tinnitus are just about
% guaranteed to hear the Zwicker tone.
%
% Lucas C. Parra, Barak A. Pearlmutter, "Illusory percepts from auditory
% adaptation", Journal of the Acoustical Society of America, accepted for
% publication November 20, 2006.
% http://newton.bme.columbia.edu/~lparra/publish/tinnitus.pdf
% 
% On a technical note. Most PC sound-cards these days, and probably a fair
% number of PC speakers will not reproduce the sharp band-gap edges nor the
% deep noise floor required to generate the Zwicker precept. We used a
% MAYA44 external USB device and the windows matlab driver for AISO
% multi-channel i/o: http://sourceforge.net/projects/pa-wavplay/ I can't
% guarentee anything else. If you use that driver you can use this code plus
% a microphone to check what you are actually putting out on your
% speakers/headphones.
%
% (c) Lucas C. Parra, November 28, 2006

if 1
    f=[0 500 1000 2000];
    bw = 2*f; % Q=2
    bw=[5000 5000 5000 5000]; 
    f0=[0 0 0 0];
else
    % lower frequency of the bandgap. use 0 for no gap
    f = [0 0 0 0 500 500 500 500 500 1000 1000 1000 1000 1000 2000 2000 2000 2000 2000]'; 
    bw = f;
    % frequency of tone to append (=0 for none, =1 for white noise)
    f0 =[0 500 1000 2000 0 1 500 1000 2000 0 1 500 1000 2000 0 1 500 1000 2000]'; 
end
a=ones(size(f0));

fs=44100;      % sampling frequency
T = fs/25;     % 40ms

% randomize order
indx=randperm(length(f));
f = f(indx);
a = a(indx);
f0 = f0(indx);

% amplitude ramp-up in fs samples and ramp-down in T samples
ru=(1:fs)'/fs;
rd=[(T:-1:1)'/T];


x=randn(fs,1)/10;        % white noise
s=randn(fs,1)*(10^-6);   % silence

stim=[];
for i=1:length(f)
    if f(i)>0
        X=fft(x); X([f(i)+(1:bw(i)) end-f(i)+2-(1:bw(i))])=0; xn=real(ifft(X));
        xn = xn*std(x)/std(xn);    
    else
        xn=x;
    end 
    switch f0(i)   
        case 0  
            stim = [stim; ru.*xn;xn;xn;xn(1:T).*rd;s;s];
        case 1     
            r = std(xn)/std(x);
            stim = [stim; ru.*xn;xn;xn;x(1:T).*rd*r;s;s];
        otherwise
            t=sin(2*pi*f0(i)*[1:fs]'/fs);
            r=std(xn)/std(t)/2;
            stim = [stim; ru.*xn;xn;xn;xn(1:T).*rd;t(1:2*T).*[rd(end:-1:1);rd]*r;s;s];
    end
end

return

if exist('pa_wavplayrecord.m')
  y = pa_wavplayrecord([stim stim],0,fs,length(stim),2,2,0,'asio'); 
else
  wavplay([stim stim],fs);
  y=stim;
end

subplot(3,1,1); specgram(stim,512,fs); title('Stimulus')
subplot(3,1,2); specgram(y,512,fs); title('Actual recorded sound')
subplot(3,1,3); plot(y); title('Actual recorded sound')
f