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
tomwalters@0: % function to generate an artificial vowel
tomwalters@0: function sig=gen_frog(sig,options)
tomwalters@0: % vowel,pitch,scale,halflifeconstant,onsettimeconstant,pitch_jitter,formant_jitter,formant_jitter_bw);% 
tomwalters@0: 
tomwalters@0: if nargin < 2
tomwalters@0:     options=[];
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % the rise of a damped sinusoid is not instantaneous, but like a gamma
tomwalters@0: % functin with the following power of t:
tomwalters@0: if isfield(options,'onsettimeconstant')
tomwalters@0:     onsettimeconstant=options.onsettimeconstant;
tomwalters@0: else
tomwalters@0:     onsettimeconstant=0.5;
tomwalters@0: end% halflife of each formant is calculated by dividing this number by the
tomwalters@0: % formant frequency:
tomwalters@0: % if halflifeconstant<=0, then fixed to 4 ms
tomwalters@0: if isfield(options,'halflifeconstant')
tomwalters@0:     halflifeconstant=options.halflifeconstant;
tomwalters@0: else
tomwalters@0:     halflifeconstant=3;
tomwalters@0: end
tomwalters@0: % scaling = 1: normal speaker
tomwalters@0: if isfield(options,'scale')
tomwalters@0:     scale=options.scale;
tomwalters@0: else
tomwalters@0:     scale=1;
tomwalters@0: end
tomwalters@0: % f0
tomwalters@0: if isfield(options,'pitch')
tomwalters@0:     pitch=options.pitch;
tomwalters@0: else
tomwalters@0:     pitch=100;
tomwalters@0: end
tomwalters@0: % how long the decay should continue longer then the reprate
tomwalters@0: if isfield(options,'carry_decay')
tomwalters@0:     carry_decay_parameter=options.carry_decay;
tomwalters@0: else
tomwalters@0:     carry_decay_parameter=1;
tomwalters@0: end
tomwalters@0: % normal or is it the octave? In this case, we jump over each second 
tomwalters@0: if isfield(options,'do_octave')
tomwalters@0:     do_octave=options.do_octave;
tomwalters@0: else
tomwalters@0:     do_octave=0;
tomwalters@0: end
tomwalters@0: % type of vowel: 'a','e','i','o','u'
tomwalters@0: if isfield(options,'vowel')
tomwalters@0:     vowel=options.vowel;
tomwalters@0: else
tomwalters@0:     vowel='a';
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if nargin < 1
tomwalters@0:     sig=signal(0.5,16000);
tomwalters@0: end
tomwalters@0: sr=getsr(sig);
tomwalters@0: dur_time=getlength(sig);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: nr_formants=4;
tomwalters@0: 
tomwalters@0: for i=1:nr_formants
tomwalters@0:     formfre(i)=pitch*i;
tomwalters@0: end
tomwalters@0: level(1)=0;
tomwalters@0: % level(2)=-18;
tomwalters@0: % level(3)=-30;
tomwalters@0: % level(4)=-38;
tomwalters@0: 
tomwalters@0: level(2)=-4;
tomwalters@0: level(3)=-6;
tomwalters@0: level(4)=-9;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % % fix all formants to the nearest harmonic of the fundamental
tomwalters@0: % if isfield(options,'adjust_to_nearest_harmonic')
tomwalters@0: %     if options.adjust_to_nearest_harmonic==1
tomwalters@0: %         f0=options.pitch;
tomwalters@0: %         for formant=1:nr_formants
tomwalters@0: %             fre=formfre(formant);
tomwalters@0: %             newfre=round(fre/f0)*f0;
tomwalters@0: %             formfre(formant)=newfre;
tomwalters@0: %         end
tomwalters@0: %     end
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: formfre=formfre.*scale;
tomwalters@0: dur_samp=dur_time.*sr;        %length of vowel defined in sample points
tomwalters@0: sgpp=1/pitch; %standard glottal pulse period.
tomwalters@0: pitch_jitter=0;
tomwalters@0: t=[0:1/sr:(dur_time-(1/sr))]; %our time sequence
tomwalters@0: dur_samp=dur_time.*sr;        %length of vowel defined in sample points
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %we now generate a sequence, specified in sample points, which
tomwalters@0: %define the spacing of glottal pulses.  With zero pitch_jitter the sequence
tomwalters@0: %is regular (1/pitch). With a pitch_jitter value of 1 the spacing of glottal 
tomwalters@0: %pulses fall in a range with an upper limit of double the period of
tomwalters@0: %the pitch and a lower limit of 1/sampling rate.
tomwalters@0: %The average spacing of glottal pulses is approx 1/pitch
tomwalters@0: lwr_pth_jit=0.5-(pitch_jitter/2);
tomwalters@0: upr_pth_jit=0.5+(pitch_jitter/2);
tomwalters@0: 
tomwalters@0: gp=cell(nr_formants,1); %for each formant will store pulse time
tomwalters@0: for formant=1:nr_formants
tomwalters@0:     enough_pulses=0;
tomwalters@0:     pulse_number=0;
tomwalters@0:     while enough_pulses==0
tomwalters@0:         pulse_number=pulse_number+1;
tomwalters@0:         pulse_spacing(pulse_number)=max(1,floor(sr.*(sgpp.*2.*(lwr_pth_jit+(upr_pth_jit-lwr_pth_jit)*rand(1)))));
tomwalters@0:         pulse_time(pulse_number)=sum(pulse_spacing)-pulse_spacing(1)+1;
tomwalters@0:         if pulse_time(pulse_number)>=dur_samp
tomwalters@0:             pulse_time=pulse_time(1:pulse_number-1);
tomwalters@0:             pulse_number=pulse_number-1;
tomwalters@0:             enough_pulses=1;
tomwalters@0:         else
tomwalters@0:         end
tomwalters@0:     end
tomwalters@0:     
tomwalters@0:     gp{formant}=pulse_time;
tomwalters@0:     clear pulse_time;
tomwalters@0:     clear pulse_spacing;
tomwalters@0:     no_pulses(formant)=length(gp{formant});
tomwalters@0:     pulse_times{formant}=gp{formant};
tomwalters@0: end%formant
tomwalters@0: 
tomwalters@0: formant_jitter=0;
tomwalters@0: 
tomwalters@0: % if no gamma tone, then precalculate the damping function once for all
tomwalters@0: % if onsettimeconstant <= 0 && halflifeconstant==0
tomwalters@0:     hl=0.04;
tomwalters@0:     damping=exp(t.*log(0.5)/hl);  %this decays to 0.5 after hl seconds
tomwalters@0:     damping=damping/max(damping);
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: onsettimes=power(t,onsettimeconstant);
tomwalters@0: 
tomwalters@0: formant_jitter_bw=0;
tomwalters@0: %--------------------------------------------------------------------------
tomwalters@0: for formant=1:nr_formants
tomwalters@0:     lw_log(formant)=max(log10(150), (log10(formfre(formant))- ((formant_jitter_bw.*.3)/2) ) ); %cannot go below 150Hz
tomwalters@0:     lw_log_adj(formant)=log10(formfre(formant))-((log10(formfre(formant))-lw_log(formant)).*formant_jitter);
tomwalters@0:     up_log(formant)=min(log10(4500),log10(formfre(formant))+ ((formant_jitter_bw.*.3)/2))  ; 
tomwalters@0:     up_log_adj(formant)=log10(formfre(formant))+((up_log(formant)-log10(formfre(formant))).*formant_jitter);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % if we want to generate the octave, we take exactly the same pulses and
tomwalters@0: % frequencies, but only every second one:
tomwalters@0: if do_octave
tomwalters@0:     pulse_step=2;
tomwalters@0: else
tomwalters@0:     pulse_step=1;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: nr_points=length(t);
tomwalters@0: final_wave=zeros(dur_samp,nr_formants);
tomwalters@0: for formant=1:nr_formants
tomwalters@0:     for count=1:pulse_step:no_pulses(formant)-1
tomwalters@0:         oneortwo=randperm(2);
tomwalters@0:         if oneortwo(1)==1 
tomwalters@0:             freq=10.^((lw_log_adj(formant)+(log10(formfre(formant))-lw_log_adj(formant))*rand(1))); %lower range
tomwalters@0:         else
tomwalters@0:             freq=10.^((log10(formfre(formant))+(up_log_adj(formant)-log10(formfre(formant)))*rand(1))); %upper range
tomwalters@0:         end
tomwalters@0: 
tomwalters@0: %         if freq>1000
tomwalters@0: %             no_octave=((log10(freq)-3))/.3;
tomwalters@0: %             lev=-12.*no_octave;
tomwalters@0: %         else
tomwalters@0: %             lev=0;
tomwalters@0: %         end
tomwalters@0: %         amp=10.^(lev/20);
tomwalters@0: 
tomwalters@0:         lev=level(formant);
tomwalters@0:         amp=10.^(lev/20);
tomwalters@0: 
tomwalters@0:         if halflifeconstant<=0
tomwalters@0:             hl=0.004;
tomwalters@0:         else
tomwalters@0:             hl=halflifeconstant/freq;
tomwalters@0:         end
tomwalters@0: 
tomwalters@0:         nr_relevant_points=pulse_times{formant}(count+1)-pulse_times{formant}(count);
tomwalters@0:         
tomwalters@0:         nr_relevant_points=nr_relevant_points*carry_decay_parameter;
tomwalters@0:         
tomwalters@0:         if onsettimeconstant > 0
tomwalters@0:             damping=zeros(nr_relevant_points,1);
tomwalters@0:             for jj=1:nr_relevant_points
tomwalters@0:                 damping(jj)=onsettimes(jj)*exp(t(jj)*log(0.5)/hl);  %this decays to 0.5 after hl seconds
tomwalters@0:             end
tomwalters@0:             damping=damping/max(damping);
tomwalters@0:         end
tomwalters@0:         grafix=0;
tomwalters@0:         if grafix==1
tomwalters@0:             figure(234)
tomwalters@0:             plot(damping);
tomwalters@0:         end
tomwalters@0:         
tomwalters@0:         this_formant=zeros(nr_relevant_points,1);
tomwalters@0:         for jj=1:nr_relevant_points
tomwalters@0:             sinsin=sin(2*pi*freq*t(jj));
tomwalters@0:             this_formant(jj)=amp*sinsin*damping(jj);
tomwalters@0:         end
tomwalters@0:         start_nr=pulse_times{formant}(count);
tomwalters@0:         stop_nr=start_nr+nr_relevant_points-1;
tomwalters@0:         
tomwalters@0:         if stop_nr> nr_points
tomwalters@0:             nr_new=nr_points-start_nr+1;
tomwalters@0:             this_formant=this_formant(1:nr_new);
tomwalters@0:             stop_nr=nr_points;
tomwalters@0:             final_wave(start_nr:stop_nr,formant)= final_wave(start_nr:stop_nr,formant)+this_formant;
tomwalters@0:         else
tomwalters@0:             final_wave(start_nr:stop_nr,formant)= final_wave(start_nr:stop_nr,formant)+this_formant;
tomwalters@0:         end
tomwalters@0:     end
tomwalters@0: end
tomwalters@0: final_wave_total=zeros(dur_samp,1);
tomwalters@0: for formant=1:nr_formants
tomwalters@0:     final_wave_total=final_wave_total+final_wave(:,formant);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % sig=signal(final_wave_total,sr);
tomwalters@0: 
tomwalters@0: if isfield(options,'rise_time')
tomwalters@0:     rise_time=options.rise_time;
tomwalters@0: else
tomwalters@0:     rise_time=0.015;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if isfield(options,'fall_time')
tomwalters@0:     fall_time=options.fall_time;
tomwalters@0: else
tomwalters@0:     fall_time=0.1;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % sigvals=
tomwalters@0: final_wave_total=gate_on_off(rise_time,fall_time,final_wave_total,sr);
tomwalters@0: sig=signal(final_wave_total,sr);
tomwalters@0: % 
tomwalters@0: % hold_time=length(sig)-rise_time-fall_time;
tomwalters@0: % attack=linspace(0,1,round(rise_time*sr));
tomwalters@0: % 
tomwalters@0: % hold=ones(1,round(hold_time*sr));
tomwalters@0: % release=linspace(1,0,round(fall_time*sr));
tomwalters@0: % envelope=[attack hold release];
tomwalters@0: % envelope=envelope(1:getnrpoints(sig));
tomwalters@0: % envelope=signal(envelope,sr);
tomwalters@0: % sig=sig*envelope;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: function amp=calc_level(freq);
tomwalters@0: %determines the level of a signal at a given frequency after it 
tomwalters@0: %has been filtered with a low-pass filter of 12dB/octave at 1kHz
tomwalters@0: %first calculate if frequency is above 1kHz and if so by how many octaves
tomwalters@0: %then multiply by slope.
tomwalters@0: if freq>1000
tomwalters@0:     no_octave=((log10(freq)-3))/.3;
tomwalters@0:     lev=-12.*no_octave;
tomwalters@0: else
tomwalters@0:     lev=0;
tomwalters@0: end
tomwalters@0: amp=10.^(lev/20);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: function output=gate_on_off(onset,offset,input,sr)
tomwalters@0: %check duration of signal is large enough for both the onset and offset
tomwalters@0: %values. 
tomwalters@0: 
tomwalters@0: sig_length_samp=length(input);
tomwalters@0: onset_length_samp=floor(onset.*sr);
tomwalters@0: offset_length_samp=floor(offset.*sr);
tomwalters@0: 
tomwalters@0: if (onset_length_samp+offset_length_samp)>sig_length_samp
tomwalters@0:     output=0;
tomwalters@0:     disp('---ERROR---  onset and offset duration too large for signal');
tomwalters@0:     return
tomwalters@0: else end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %generate onset and offset amplitudes
tomwalters@0: n_on=onset_length_samp;
tomwalters@0: k=[1:1:n_on];
tomwalters@0: onset_win=(1-cos(2.*pi.*(k./(2.*(n_on-1)))))./2;
tomwalters@0: 
tomwalters@0: n_off=offset_length_samp;
tomwalters@0: k=[1:1:n_off];
tomwalters@0: offset_win=0.5+((cos(2.*pi.*(k./(2.*(n_off-1)))))./2);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: total_window=ones(sig_length_samp,1);
tomwalters@0: total_window(1:onset_length_samp)=onset_win;
tomwalters@0: total_window(sig_length_samp-offset_length_samp+1:sig_length_samp)=offset_win;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: output=input.*total_window;
tomwalters@0: 
tomwalters@0: return