bleeck@3: % method of class @signal
bleeck@3: % function sig=generatedampsinus(sig,carfre,modfre,amplitude,halflife)
bleeck@3: %   INPUT VALUES:
bleeck@3: %       sig: original @signal with length and samplerate
bleeck@3: %       carfre: carrier frequency (Hz) [1000]
bleeck@3: %       modfre: modulation frequency (Hz) [100]
bleeck@3: %       amplitude: [1]
bleeck@3: %       halflife: time for the envelope envelope to decrease exponentielly
bleeck@3: %       to 1/2
bleeck@3: %
bleeck@3: %
bleeck@3: %   RETURN VALUE:
bleeck@3: %       sig:  @signal
bleeck@3: %
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: 
bleeck@3: 
bleeck@3: 
bleeck@3: function sig=generatedampsinus(sig,carfre,modfre,amplitude,halflife,jitter)
bleeck@3: % generates a damped sinusoid, that is a carrier pure tone modulated with a
bleeck@3: % exponentially decreasing envelope.
bleeck@3: % sig is the signal
bleeck@3: % carfre is the carrier frequency of the pure tone
bleeck@3: % modfre is the modulation frequency (in Hz)
bleeck@3: % amplitude is the final amplitude
bleeck@3: % halflife is the time in seconds, in which the envelope drops to its half value
bleeck@3: % if halflife is 'gamma', instead a gammaenvelope is used
bleeck@3: 
bleeck@3: % jitter is how regular the pulses are (0-1) 1=100%
bleeck@3: % if jitter is 'repeat', only one pulse is generated and repeated
bleeck@3: 
bleeck@3: 
bleeck@3: if nargin < 6
bleeck@3:     jitter=0;
bleeck@3: end
bleeck@3: if nargin < 5
bleeck@3:     halflife=0.01;
bleeck@3: end
bleeck@3: if nargin < 4
bleeck@3:     amplitude=1;
bleeck@3: end
bleeck@3: 
bleeck@3: if nargin < 3
bleeck@3:     modfre=100;
bleeck@3: end
bleeck@3: if nargin < 2
bleeck@3:     carfre=1000;
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: if isequal(halflife,'gamma') &&  isequal(jitter,'repeat')
bleeck@3:     % code for Tom's 2010 MSc project
bleeck@3:     
bleeck@3:     % generate only one period and repeat it n times
bleeck@3:     sr=getsr(sig);
bleeck@3:     periodlen=1/modfre;
bleeck@3:     oneperiod=signal(periodlen,sr);
bleeck@3:     
bleeck@3:     sinus=generatesinus(oneperiod,carfre,amplitude,0);
bleeck@3:     
bleeck@3:     % calculate envelope and mult both
bleeck@3:     envelope=oneperiod;
bleeck@3:     t=[0:1/sr:periodlen-1/sr]*500*(500/getnrpoints(envelope));
bleeck@3:     env=power(t,4).*exp(-2*pi*t/2);
bleeck@3:     env=env./max(env);
bleeck@3:     envelope=setvalues(envelope,env);
bleeck@3:     envelope=envelope/max(envelope)*amplitude;
bleeck@3:     envelope=setstarttime(envelope,0);
bleeck@3:     
bleeck@3:     % set the envelope and the amplitude
bleeck@3:     oneperiod=sinus*envelope;
bleeck@3:     
bleeck@3:     % repeat for the right number of repeats
bleeck@3:     nr_repeats=round(getlength(sig)/getlength(envelope));
bleeck@3:     fsig=oneperiod;
bleeck@3:     
bleeck@3:     if nr_repeats>1
bleeck@3:         for i=1:nr_repeats-1
bleeck@3:             fsig=append(fsig,oneperiod);
bleeck@3:         end
bleeck@3:     end
bleeck@3:     
bleeck@3:     sig=fsig;
bleeck@3:     
bleeck@3:     sig=setname(sig,sprintf('damped sinusoid %4.2f kHz, Modulation=%4.1f Hz, gamma envelope',carfre/1000,modfre));
bleeck@3:     
bleeck@3:     %     figure(1)
bleeck@3:     %     plot(env)
bleeck@3:     
bleeck@3:     return
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: sinus=generatesinus(sig,carfre,amplitude,0);
bleeck@3: 
bleeck@3: % calculate envelope and mult both
bleeck@3: envelope=sig;
bleeck@3: 
bleeck@3: time_const=halflife/0.69314718;
bleeck@3: 
bleeck@3: env_vals=getvalues(envelope);
bleeck@3: sr=getsr(envelope);
bleeck@3: reprate=1/modfre;
bleeck@3: 
bleeck@3: 
bleeck@3: sig_len=getlength(sig);
bleeck@3: % when regular, all modulataions are at the same time:
bleeck@3: if jitter==0
bleeck@3:     pulse_times=0:reprate:sig_len;
bleeck@3: else
bleeck@3:     nr_pulses=ceil(sig_len/reprate);
bleeck@3:     modulation_period=1/modfre;
bleeck@3:     pulse_times(1)=0;
bleeck@3:     for n = 2:nr_pulses
bleeck@3:         jittering=(rand-0.5*jitter)*modulation_period;
bleeck@3:         pulse_times(n) = pulse_times(n-1) + modulation_period+jittering;
bleeck@3:     end
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: % oldval=0;
bleeck@3: % corr=exp(1/sr/time_const);
bleeck@3: % t=0;
bleeck@3: next_pulse=pulse_times(1);
bleeck@3: pulse_counter=1;
bleeck@3: 
bleeck@3: % for i=1:getnrpoints(envelope);
bleeck@3: %     oldval=oldval/corr;
bleeck@3: %     t=t+1/sr;
bleeck@3: %     if t>next_pulse
bleeck@3: %         oldval=1;
bleeck@3: %         pulse_counter=pulse_counter+1;
bleeck@3: %         if length(pulse_times)>=pulse_counter
bleeck@3: %             next_pulse=pulse_times(pulse_counter);
bleeck@3: %         else
bleeck@3: %             next_pulse=inf;
bleeck@3: %         end
bleeck@3: %     end
bleeck@3: %     env_vals(i)= oldval;
bleeck@3: % end
bleeck@3: 
bleeck@3: % onsettimeconstant
bleeck@3: %
bleeck@3: t=[0:1/sr:reprate-1/sr]*500;
bleeck@3: env=power(t,4).*exp(-2*pi*t/2);
bleeck@3: env=env./max(env);
bleeck@3: 
bleeck@3: % e=linspace(env(length(env)-l),0,l+1);
bleeck@3: % env(length(env)-l:end)=e;
bleeck@3: 
bleeck@3: %  e=linspace(1,0,l+1);
bleeck@3: % env(length(env)-l:end)=env(length(env)-l:end).*e.*e;
bleeck@3: 
bleeck@3: t=0;
bleeck@3: ct=1;
bleeck@3: for i=1:getnrpoints(envelope);
bleeck@3:     t=t+1/sr;
bleeck@3:     ct=ct+1;
bleeck@3:     if t>next_pulse
bleeck@3:         ct=1;
bleeck@3:         pulse_counter=pulse_counter+1;
bleeck@3:         if length(pulse_times)>=pulse_counter
bleeck@3:             next_pulse=pulse_times(pulse_counter);
bleeck@3:         else
bleeck@3:             next_pulse=inf;
bleeck@3:         end
bleeck@3:     end
bleeck@3:     if ct>length(env)
bleeck@3:         env_vals(i)= 0;
bleeck@3:     else
bleeck@3:         env_vals(i)= env(ct);
bleeck@3:     end
bleeck@3: end
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: %
bleeck@3: % onsettime=0;
bleeck@3: % for i=1:getnrpoints(envelope);
bleeck@3: %     t=t+1/sr;
bleeck@3: %
bleeck@3: % %     env_vals(i)= exp(-(time)/time_const);
bleeck@3: %     env_vals(i)= power(t,onsettime)*exp(-(t)/time_const);
bleeck@3: %     t=mod(t,reprate);
bleeck@3: %
bleeck@3: % end
bleeck@3: envelope=setvalues(envelope,env_vals);
bleeck@3: envelope=envelope/max(envelope)*amplitude;
bleeck@3: envelope=setstarttime(envelope,0);
bleeck@3: 
bleeck@3: % set the envelope and the amplitude
bleeck@3: sig=sinus*envelope;
bleeck@3: 
bleeck@3: sig=setname(sig,sprintf('damped sinusoid %4.2f kHz, Modulation=%4.1f Hz, halflife=%4.1f ms',carfre/1000,modfre,halflife*1000));
bleeck@3: 
bleeck@3: % plot(sig)
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: 
bleeck@3: %
bleeck@3: %
bleeck@3: %
bleeck@3: % % from here on old code, might not work!
bleeck@3: %
bleeck@3: % sinus=generatesinus(sig,carfre,amplitude,0);
bleeck@3: %
bleeck@3: % % calculate envelope and mult both
bleeck@3: % envelope=sig;
bleeck@3: %
bleeck@3: % time_const=halflife/0.69314718;
bleeck@3: %
bleeck@3: % env_vals=getvalues(envelope);
bleeck@3: % sr=getsr(envelope);
bleeck@3: % reprate=1/modfre;
bleeck@3: %
bleeck@3: %
bleeck@3: % sig_len=getlength(sig);
bleeck@3: % % when regular, all modulataions are at the same time:
bleeck@3: % if jitter==0
bleeck@3: %     pulse_times=0:reprate:sig_len;
bleeck@3: % else
bleeck@3: %     nr_pulses=ceil(sig_len/reprate);
bleeck@3: %     modulation_period=1/modfre;
bleeck@3: %     pulse_times(1)=0;
bleeck@3: %     for n = 2:nr_pulses
bleeck@3: %         jittering=(rand-0.5*jitter)*modulation_period;
bleeck@3: %         pulse_times(n) = pulse_times(n-1) + modulation_period+jittering;
bleeck@3: %     end
bleeck@3: % end
bleeck@3: %
bleeck@3: %
bleeck@3: %
bleeck@3: % % oldval=0;
bleeck@3: % % corr=exp(1/sr/time_const);
bleeck@3: % % t=0;
bleeck@3: % next_pulse=pulse_times(1);
bleeck@3: % pulse_counter=1;
bleeck@3: %
bleeck@3: % % for i=1:getnrpoints(envelope);
bleeck@3: % %     oldval=oldval/corr;
bleeck@3: % %     t=t+1/sr;
bleeck@3: % %     if t>next_pulse
bleeck@3: % %         oldval=1;
bleeck@3: % %         pulse_counter=pulse_counter+1;
bleeck@3: % %         if length(pulse_times)>=pulse_counter
bleeck@3: % %             next_pulse=pulse_times(pulse_counter);
bleeck@3: % %         else
bleeck@3: % %             next_pulse=inf;
bleeck@3: % %         end
bleeck@3: % %     end
bleeck@3: % %     env_vals(i)= oldval;
bleeck@3: % % end
bleeck@3: %
bleeck@3: % % onsettimeconstant
bleeck@3: % %
bleeck@3: % t=[0:1/sr:reprate-1/sr]*500;
bleeck@3: % env=power(t,4).*exp(-2*pi*t/2);
bleeck@3: % env=env./max(env);
bleeck@3: %
bleeck@3: % % e=linspace(env(length(env)-l),0,l+1);
bleeck@3: % % env(length(env)-l:end)=e;
bleeck@3: %
bleeck@3: % %  e=linspace(1,0,l+1);
bleeck@3: % % env(length(env)-l:end)=env(length(env)-l:end).*e.*e;
bleeck@3: %
bleeck@3: % t=0;
bleeck@3: % ct=1;
bleeck@3: % for i=1:getnrpoints(envelope);
bleeck@3: %     t=t+1/sr;
bleeck@3: %     ct=ct+1;
bleeck@3: %     if t>next_pulse
bleeck@3: %         ct=1;
bleeck@3: %         pulse_counter=pulse_counter+1;
bleeck@3: %         if length(pulse_times)>=pulse_counter
bleeck@3: %             next_pulse=pulse_times(pulse_counter);
bleeck@3: %         else
bleeck@3: %             next_pulse=inf;
bleeck@3: %         end
bleeck@3: %     end
bleeck@3: %     if ct>length(env)
bleeck@3: %         env_vals(i)= 0;
bleeck@3: %     else
bleeck@3: %         env_vals(i)= env(ct);
bleeck@3: %     end
bleeck@3: % end
bleeck@3: %
bleeck@3: %
bleeck@3: %
bleeck@3: % %
bleeck@3: % % onsettime=0;
bleeck@3: % % for i=1:getnrpoints(envelope);
bleeck@3: % %     t=t+1/sr;
bleeck@3: % %
bleeck@3: % % %     env_vals(i)= exp(-(time)/time_const);
bleeck@3: % %     env_vals(i)= power(t,onsettime)*exp(-(t)/time_const);
bleeck@3: % %     t=mod(t,reprate);
bleeck@3: % %
bleeck@3: % % end
bleeck@3: % envelope=setvalues(envelope,env_vals);
bleeck@3: % envelope=envelope/max(envelope)*amplitude;
bleeck@3: % envelope=setstarttime(envelope,0);
bleeck@3: %
bleeck@3: % % set the envelope and the amplitude
bleeck@3: % sig=sinus*envelope;
bleeck@3: %
bleeck@3: % sig=setname(sig,sprintf('damped sinusoid %4.2f kHz, Modulation=%4.1f Hz, halflife=%4.1f ms',carfre/1000,modfre,halflife*1000));
bleeck@3: %
bleeck@3: % % plot(sig)