bleeck@3: % This external file is included as part of the 'aim-mat' distribution package
bleeck@3: % (c) 2011, University of Southampton
bleeck@3: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@3: % download of current version is on the soundsoftware site: 
bleeck@3: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@3: % documentation and everything is on http://www.acousticscale.org
bleeck@3: % function to generate an artificial vowel
bleeck@3: function sig=gen_vowel(sig,options)
bleeck@3: % vowel,pitch,scale,halflifeconstant,onsettimeconstant,pitch_jitter,formant_jitter,formant_jitter_bw);% 
bleeck@3: 
bleeck@3: if nargin < 2
bleeck@3:     options=[];
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: % jitterwidth of jitter in formant frequencies
bleeck@3: %formant_jitter_bw=1; %bandwidth specified in octaves over which the formants are jittered
bleeck@3: if isfield(options,'formant_jitter_bw')
bleeck@3:     formant_jitter_bw=options.formant_jitter_bw;
bleeck@3: else
bleeck@3:     formant_jitter_bw=1; % in octaves
bleeck@3: end
bleeck@3: % jitter in percent of bandwidth
bleeck@3: %formant_jitter=0.0; %------- range is 0 to 1, representing 0% to 100%
bleeck@3: if isfield(options,'formant_jitter')
bleeck@3:     formant_jitter=options.formant_jitter;
bleeck@3: else
bleeck@3:     formant_jitter=0.03; % 
bleeck@3: end
bleeck@3: % jitter in pitch
bleeck@3: %pitch_jitter=0.0; %------- range is 0 to 1, representing 0% to 100%
bleeck@3: if isfield(options,'pitch_jitter')
bleeck@3:     pitch_jitter=options.pitch_jitter;
bleeck@3: else
bleeck@3:     pitch_jitter=0.01;
bleeck@3: end
bleeck@3: % the rise of a damped sinusoid is not instantaneous, but like a gamma
bleeck@3: % functin with the following power of t:
bleeck@3: if isfield(options,'onsettimeconstant')
bleeck@3:     onsettimeconstant=options.onsettimeconstant;
bleeck@3: else
bleeck@3:     onsettimeconstant=0.5;
bleeck@3: end% halflife of each formant is calculated by dividing this number by the
bleeck@3: % formant frequency:
bleeck@3: % if halflifeconstant<=0, then fixed to 4 ms
bleeck@3: if isfield(options,'halflifeconstant')
bleeck@3:     halflifeconstant=options.halflifeconstant;
bleeck@3: else
bleeck@3:     halflifeconstant=3;
bleeck@3: end
bleeck@3: % scaling = 1: normal speaker
bleeck@3: if isfield(options,'scale')
bleeck@3:     scale=options.scale;
bleeck@3: else
bleeck@3:     scale=1;
bleeck@3: end
bleeck@3: % f0
bleeck@3: if isfield(options,'pitch')
bleeck@3:     pitch=options.pitch;
bleeck@3: else
bleeck@3:     pitch=100;
bleeck@3: end
bleeck@3: % how long the decay should continue longer then the reprate
bleeck@3: if isfield(options,'carry_decay')
bleeck@3:     carry_decay_parameter=options.carry_decay;
bleeck@3: else
bleeck@3:     carry_decay_parameter=1;
bleeck@3: end
bleeck@3: % normal or is it the octave? In this case, we jump over each second 
bleeck@3: if isfield(options,'do_octave')
bleeck@3:     do_octave=options.do_octave;
bleeck@3: else
bleeck@3:     do_octave=0;
bleeck@3: end
bleeck@3: % type of vowel: 'a','e','i','o','u'
bleeck@3: if isfield(options,'vowel')
bleeck@3:     vowel=options.vowel;
bleeck@3: else
bleeck@3:     vowel='a';
bleeck@3: end
bleeck@3: 
bleeck@3: % 
bleeck@3: % which formants
bleeck@3: if isfield(options,'use_formants')
bleeck@3:     use_formants=options.use_formants;
bleeck@3: else
bleeck@3:     use_formants=[1 2 3 4];
bleeck@3: end
bleeck@3: 
bleeck@3: if nargin < 1
bleeck@3:     sig=signal(0.5,16000);
bleeck@3: end
bleeck@3: sr=getsr(sig);
bleeck@3: dur_time=getlength(sig);
bleeck@3: 
bleeck@3: 
bleeck@3: nr_formants=length(use_formants);
bleeck@3: 
bleeck@3: 
bleeck@3: switch vowel
bleeck@3:     case 'a';
bleeck@3:         formfre_org(1)=730; formfre_org(2)=1090; formfre_org(3)=2440; formfre_org(4)=3300;% standard formant freqs       
bleeck@3:     case 'e';
bleeck@3:         formfre_org(1)=400; formfre_org(2)=1990; formfre_org(3)=2550; formfre_org(4)=3700; % standard formant freqs       
bleeck@3:     case 'i';
bleeck@3:         formfre_org(1)=270; formfre_org(2)=2260; formfre_org(3)=3000; formfre_org(4)=3700;% standard formant freqs     
bleeck@3:     case 'o';
bleeck@3:         formfre_org(1)=570; formfre_org(2)=850; formfre_org(3)=2430; formfre_org(4)=3300;% standard formant freqs
bleeck@3:     case 'u';
bleeck@3:         formfre_org(1)=300; formfre_org(2)=850; formfre_org(3)=2260; formfre_org(4)=3700;% standard formant freqs   
bleeck@3:     otherwise
bleeck@3:         error('signal::gen_vowel: dont know vowel');
bleeck@3: end
bleeck@3: 
bleeck@3: for i=1:length(use_formants)
bleeck@3:     formfre(i)=formfre_org(use_formants(i));
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: % fix all formants to the nearest harmonic of the fundamental
bleeck@3: if isfield(options,'adjust_to_nearest_harmonic')
bleeck@3:     if options.adjust_to_nearest_harmonic==1
bleeck@3:         f0=options.pitch;
bleeck@3:         for formant=1:nr_formants
bleeck@3:             fre=formfre(formant);
bleeck@3:             newfre=round(fre/f0)*f0;
bleeck@3:             formfre(formant)=newfre;
bleeck@3:         end
bleeck@3:     end
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: formfre=formfre.*scale;
bleeck@3: dur_samp=dur_time.*sr;        %length of vowel defined in sample points
bleeck@3: sgpp=1/pitch; %standard glottal pulse period.
bleeck@3: pitch_jitter=max(0,min(1,pitch_jitter));
bleeck@3: t=[0:1/sr:(dur_time-(1/sr))]; %our time sequence
bleeck@3: dur_samp=dur_time.*sr;        %length of vowel defined in sample points
bleeck@3: 
bleeck@3: 
bleeck@3: %we now generate a sequence, specified in sample points, which
bleeck@3: %define the spacing of glottal pulses.  With zero pitch_jitter the sequence
bleeck@3: %is regular (1/pitch). With a pitch_jitter value of 1 the spacing of glottal 
bleeck@3: %pulses fall in a range with an upper limit of double the period of
bleeck@3: %the pitch and a lower limit of 1/sampling rate.
bleeck@3: %The average spacing of glottal pulses is approx 1/pitch
bleeck@3: lwr_pth_jit=0.5-(pitch_jitter/2);
bleeck@3: upr_pth_jit=0.5+(pitch_jitter/2);
bleeck@3: 
bleeck@3: gp=cell(nr_formants,1); %for each formant will store pulse time
bleeck@3: for formant=1:nr_formants
bleeck@3:     enough_pulses=0;
bleeck@3:     pulse_number=0;
bleeck@3:     while enough_pulses==0
bleeck@3:         pulse_number=pulse_number+1;
bleeck@3:         pulse_spacing(pulse_number)=max(1,floor(sr.*(sgpp.*2.*(lwr_pth_jit+(upr_pth_jit-lwr_pth_jit)*rand(1)))));
bleeck@3:         pulse_time(pulse_number)=sum(pulse_spacing)-pulse_spacing(1)+1;
bleeck@3:         if pulse_time(pulse_number)>=dur_samp
bleeck@3:             pulse_time=pulse_time(1:pulse_number-1);
bleeck@3:             pulse_number=pulse_number-1;
bleeck@3:             enough_pulses=1;
bleeck@3:         else
bleeck@3:         end
bleeck@3:     end
bleeck@3:     
bleeck@3:     
bleeck@3:     gp{formant}=pulse_time;
bleeck@3:     clear pulse_time;
bleeck@3:     clear pulse_spacing;
bleeck@3:     no_pulses(formant)=length(gp{formant});
bleeck@3:     pulse_times{formant}=gp{formant};
bleeck@3: end%formant
bleeck@3: 
bleeck@3: formant_jitter=max(0,min(1,formant_jitter));
bleeck@3: 
bleeck@3: % if no gamma tone, then precalculate the damping function once for all
bleeck@3: % if onsettimeconstant <= 0 && halflifeconstant==0
bleeck@3:     hl=0.04;
bleeck@3:     damping=exp(t.*log(0.5)/hl);  %this decays to 0.5 after hl seconds
bleeck@3:     damping=damping/max(damping);
bleeck@3: % end
bleeck@3: 
bleeck@3: onsettimes=power(t,onsettimeconstant);
bleeck@3: 
bleeck@3: 
bleeck@3: %--------------------------------------------------------------------------
bleeck@3: for formant=1:nr_formants
bleeck@3:     lw_log(formant)=max(log10(150), (log10(formfre(formant))- ((formant_jitter_bw.*.3)/2) ) ); %cannot go below 150Hz
bleeck@3:     lw_log_adj(formant)=log10(formfre(formant))-((log10(formfre(formant))-lw_log(formant)).*formant_jitter);
bleeck@3:     up_log(formant)=min(log10(4500),log10(formfre(formant))+ ((formant_jitter_bw.*.3)/2))  ; 
bleeck@3:     up_log_adj(formant)=log10(formfre(formant))+((up_log(formant)-log10(formfre(formant))).*formant_jitter);
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: % if we want to generate the octave, we take exactly the same pulses and
bleeck@3: % frequencies, but only every second one:
bleeck@3: if do_octave
bleeck@3:     pulse_step=2;
bleeck@3: else
bleeck@3:     pulse_step=1;
bleeck@3: end
bleeck@3: 
bleeck@3: nr_points=length(t);
bleeck@3: final_wave=zeros(dur_samp,nr_formants);
bleeck@3: for formant=1:nr_formants
bleeck@3:     for count=1:pulse_step:no_pulses(formant)-1
bleeck@3:         oneortwo=randperm(2);
bleeck@3:         if oneortwo(1)==1 
bleeck@3:             freq=10.^((lw_log_adj(formant)+(log10(formfre(formant))-lw_log_adj(formant))*rand(1))); %lower range
bleeck@3:         else
bleeck@3:             freq=10.^((log10(formfre(formant))+(up_log_adj(formant)-log10(formfre(formant)))*rand(1))); %upper range
bleeck@3:         end
bleeck@3: 
bleeck@3:         if freq>1000
bleeck@3:             no_octave=((log10(freq)-3))/.3;
bleeck@3:             lev=-12.*no_octave;
bleeck@3:         else
bleeck@3:             lev=0;
bleeck@3:         end
bleeck@3:         amp=10.^(lev/20);
bleeck@3:         
bleeck@3:         if halflifeconstant<=0
bleeck@3:             hl=0.004;
bleeck@3:         else
bleeck@3:             hl=halflifeconstant/freq;
bleeck@3:         end
bleeck@3: 
bleeck@3:         nr_relevant_points=pulse_times{formant}(count+1)-pulse_times{formant}(count);
bleeck@3:         
bleeck@3:         nr_relevant_points=nr_relevant_points*carry_decay_parameter;
bleeck@3:         
bleeck@3:         if onsettimeconstant > 0
bleeck@3:             damping=zeros(nr_relevant_points,1);
bleeck@3:             for jj=1:nr_relevant_points
bleeck@3:                 damping(jj)=onsettimes(jj)*exp(t(jj)*log(0.5)/hl);  %this decays to 0.5 after hl seconds
bleeck@3:             end
bleeck@3:             damping=damping/max(damping);
bleeck@3:         end
bleeck@3:         grafix=0;
bleeck@3:         if grafix==1
bleeck@3:             figure(234)
bleeck@3:             plot(damping);
bleeck@3:         end
bleeck@3:         
bleeck@3:         this_formant=zeros(nr_relevant_points,1);
bleeck@3:         for jj=1:nr_relevant_points
bleeck@3:             sinsin=sin(2*pi*freq*t(jj));
bleeck@3:             this_formant(jj)=amp*sinsin*damping(jj);
bleeck@3:         end
bleeck@3:         start_nr=pulse_times{formant}(count);
bleeck@3:         stop_nr=start_nr+nr_relevant_points-1;
bleeck@3:         
bleeck@3:         if stop_nr> nr_points
bleeck@3:             nr_new=nr_points-start_nr+1;
bleeck@3:             this_formant=this_formant(1:nr_new);
bleeck@3:             stop_nr=nr_points;
bleeck@3:             final_wave(start_nr:stop_nr,formant)= final_wave(start_nr:stop_nr,formant)+this_formant;
bleeck@3:         else
bleeck@3:             final_wave(start_nr:stop_nr,formant)= final_wave(start_nr:stop_nr,formant)+this_formant;
bleeck@3:         end
bleeck@3:     end
bleeck@3: end
bleeck@3: final_wave_total=zeros(dur_samp,1);
bleeck@3: for formant=1:nr_formants
bleeck@3:     final_wave_total=final_wave_total+final_wave(:,formant);
bleeck@3: end
bleeck@3: 
bleeck@3: % sig=signal(final_wave_total,sr);
bleeck@3: 
bleeck@3: if isfield(options,'rise_time')
bleeck@3:     rise_time=options.rise_time;
bleeck@3: else
bleeck@3:     rise_time=0.015;
bleeck@3: end
bleeck@3: 
bleeck@3: if isfield(options,'fall_time')
bleeck@3:     fall_time=options.fall_time;
bleeck@3: else
bleeck@3:     fall_time=0.1;
bleeck@3: end
bleeck@3: 
bleeck@3: % sigvals=
bleeck@3: final_wave_total=gate_on_off(rise_time,fall_time,final_wave_total,sr);
bleeck@3: sig=signal(final_wave_total,sr);
bleeck@3: % 
bleeck@3: % hold_time=length(sig)-rise_time-fall_time;
bleeck@3: % attack=linspace(0,1,round(rise_time*sr));
bleeck@3: % 
bleeck@3: % hold=ones(1,round(hold_time*sr));
bleeck@3: % release=linspace(1,0,round(fall_time*sr));
bleeck@3: % envelope=[attack hold release];
bleeck@3: % envelope=envelope(1:getnrpoints(sig));
bleeck@3: % envelope=signal(envelope,sr);
bleeck@3: % sig=sig*envelope;
bleeck@3: 
bleeck@3: 
bleeck@3: function amp=calc_level(freq);
bleeck@3: %determines the level of a signal at a given frequency after it 
bleeck@3: %has been filtered with a low-pass filter of 12dB/octave at 1kHz
bleeck@3: %first calculate if frequency is above 1kHz and if so by how many octaves
bleeck@3: %then multiply by slope.
bleeck@3: if freq>1000
bleeck@3:     no_octave=((log10(freq)-3))/.3;
bleeck@3:     lev=-12.*no_octave;
bleeck@3: else
bleeck@3:     lev=0;
bleeck@3: end
bleeck@3: amp=10.^(lev/20);
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: function output=gate_on_off(onset,offset,input,sr)
bleeck@3: %check duration of signal is large enough for both the onset and offset
bleeck@3: %values. 
bleeck@3: 
bleeck@3: sig_length_samp=length(input);
bleeck@3: onset_length_samp=floor(onset.*sr);
bleeck@3: offset_length_samp=floor(offset.*sr);
bleeck@3: 
bleeck@3: if (onset_length_samp+offset_length_samp)>sig_length_samp
bleeck@3:     output=0;
bleeck@3:     disp('---ERROR---  onset and offset duration too large for signal');
bleeck@3:     return
bleeck@3: else end
bleeck@3: 
bleeck@3: 
bleeck@3: %generate onset and offset amplitudes
bleeck@3: n_on=onset_length_samp;
bleeck@3: k=[1:1:n_on];
bleeck@3: onset_win=(1-cos(2.*pi.*(k./(2.*(n_on-1)))))./2;
bleeck@3: 
bleeck@3: n_off=offset_length_samp;
bleeck@3: k=[1:1:n_off];
bleeck@3: offset_win=0.5+((cos(2.*pi.*(k./(2.*(n_off-1)))))./2);
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: total_window=ones(sig_length_samp,1);
bleeck@3: total_window(1:onset_length_samp)=onset_win;
bleeck@3: total_window(sig_length_samp-offset_length_samp+1:sig_length_samp)=offset_win;
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: output=input.*total_window;
bleeck@3: 
bleeck@3: return