tomwalters@0: % method of class @signal tomwalters@0: % function sig=generatedampsinus(sig,carfre,modfre,amplitude,halflife) tomwalters@0: % INPUT VALUES: tomwalters@0: % sig: original @signal with length and samplerate tomwalters@0: % carfre: carrier frequency (Hz) [1000] tomwalters@0: % modfre: modulation frequency (Hz) [100] tomwalters@0: % amplitude: [1] tomwalters@0: % halflife: time for the envelope envelope to decrease exponentielly tomwalters@0: % to 1/2 tomwalters@0: % tomwalters@0: % tomwalters@0: % RETURN VALUE: tomwalters@0: % sig: @signal tomwalters@0: % tomwalters@0: % (c) 2003, University of Cambridge, Medical Research Council tomwalters@0: % Stefan Bleeck (stefan@bleeck.de) tomwalters@0: % http://www.mrc-cbu.cam.ac.uk/cnbh/aimmanual tomwalters@0: % $Date: 2003/07/28 14:20:20 $ tomwalters@0: % $Revision: 1.6 $ tomwalters@0: tomwalters@0: tomwalters@0: function sig=generatedampsinus(sig,carfre,modfre,amplitude,halflife,jitter) tomwalters@0: % generates a damped sinusoid, that is a carrier pure tone modulated with a tomwalters@0: % exponentially decreasing envelope. tomwalters@0: % sig is the signal tomwalters@0: % carfre is the carrier frequency of the pure tone tomwalters@0: % modfre is the modulation frequency (in Hz) tomwalters@0: % amplitude is the final amplitude tomwalters@0: % halflife is the time in seconds, in which the envelope drops to its half value tomwalters@0: % if halflife is 'gamma', instead a gammaenvelope is used tomwalters@0: tomwalters@0: % jitter is how regular the pulses are (0-1) 1=100% tomwalters@0: % if jitter is 'repeat', only one pulse is generated and repeated tomwalters@0: tomwalters@0: tomwalters@0: if nargin < 6 tomwalters@0: jitter=0; tomwalters@0: end tomwalters@0: if nargin < 5 tomwalters@0: halflife=0.01; tomwalters@0: end tomwalters@0: if nargin < 4 tomwalters@0: amplitude=1; tomwalters@0: end tomwalters@0: tomwalters@0: if nargin < 3 tomwalters@0: modfre=100; tomwalters@0: end tomwalters@0: if nargin < 2 tomwalters@0: carfre=1000; tomwalters@0: end tomwalters@0: tomwalters@0: tomwalters@0: if isequal(halflife,'gamma') && isequal(jitter,'repeat') tomwalters@0: % code for Tom's 2010 MSc project tomwalters@0: tomwalters@0: sinus=generatesinus(sig,carfre,amplitude,0); tomwalters@0: % calculate envelope and mult both tomwalters@0: envelope=sig; tomwalters@0: tomwalters@0: env_vals=getvalues(envelope); tomwalters@0: sr=getsr(envelope); tomwalters@0: reprate=1/modfre; tomwalters@0: tomwalters@0: tomwalters@0: sig_len=getlength(sig); tomwalters@0: % when regular, all modulataions are at the same time: tomwalters@0: if jitter==0 tomwalters@0: pulse_times=0:reprate:sig_len; tomwalters@0: else tomwalters@0: nr_pulses=ceil(sig_len/reprate); tomwalters@0: modulation_period=1/modfre; tomwalters@0: pulse_times(1)=0; tomwalters@0: for n = 2:nr_pulses tomwalters@0: jittering=(rand-0.5*jitter)*modulation_period; tomwalters@0: pulse_times(n) = pulse_times(n-1) + modulation_period+jittering; tomwalters@0: end tomwalters@0: end tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: % oldval=0; tomwalters@0: % corr=exp(1/sr/time_const); tomwalters@0: % t=0; tomwalters@0: next_pulse=pulse_times(1); tomwalters@0: pulse_counter=1; tomwalters@0: tomwalters@0: % for i=1:getnrpoints(envelope); tomwalters@0: % oldval=oldval/corr; tomwalters@0: % t=t+1/sr; tomwalters@0: % if t>next_pulse tomwalters@0: % oldval=1; tomwalters@0: % pulse_counter=pulse_counter+1; tomwalters@0: % if length(pulse_times)>=pulse_counter tomwalters@0: % next_pulse=pulse_times(pulse_counter); tomwalters@0: % else tomwalters@0: % next_pulse=inf; tomwalters@0: % end tomwalters@0: % end tomwalters@0: % env_vals(i)= oldval; tomwalters@0: % end tomwalters@0: tomwalters@0: % onsettimeconstant tomwalters@0: % tomwalters@0: t=[0:1/sr:reprate-1/sr]*500; tomwalters@0: env=power(t,4).*exp(-2*pi*t/2); tomwalters@0: env=env./max(env); tomwalters@0: tomwalters@0: % e=linspace(env(length(env)-l),0,l+1); tomwalters@0: % env(length(env)-l:end)=e; tomwalters@0: tomwalters@0: % e=linspace(1,0,l+1); tomwalters@0: % env(length(env)-l:end)=env(length(env)-l:end).*e.*e; tomwalters@0: tomwalters@0: t=0; tomwalters@0: ct=1; tomwalters@0: for i=1:getnrpoints(envelope); tomwalters@0: t=t+1/sr; tomwalters@0: ct=ct+1; tomwalters@0: if t>next_pulse tomwalters@0: ct=1; tomwalters@0: pulse_counter=pulse_counter+1; tomwalters@0: if length(pulse_times)>=pulse_counter tomwalters@0: next_pulse=pulse_times(pulse_counter); tomwalters@0: else tomwalters@0: next_pulse=inf; tomwalters@0: end tomwalters@0: end tomwalters@0: if ct>length(env) tomwalters@0: env_vals(i)= 0; tomwalters@0: else tomwalters@0: env_vals(i)= env(ct); tomwalters@0: end tomwalters@0: end tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: % tomwalters@0: % onsettime=0; tomwalters@0: % for i=1:getnrpoints(envelope); tomwalters@0: % t=t+1/sr; tomwalters@0: % tomwalters@0: % % env_vals(i)= exp(-(time)/time_const); tomwalters@0: % env_vals(i)= power(t,onsettime)*exp(-(t)/time_const); tomwalters@0: % t=mod(t,reprate); tomwalters@0: % tomwalters@0: % end tomwalters@0: envelope=setvalues(envelope,env_vals); tomwalters@0: envelope=envelope/max(envelope)*amplitude; tomwalters@0: envelope=setstarttime(envelope,0); tomwalters@0: tomwalters@0: % set the envelope and the amplitude tomwalters@0: sig=sinus*envelope; tomwalters@0: tomwalters@0: sig=setname(sig,sprintf('damped sinusoid %4.2f kHz, Modulation=%4.1f Hz, halflife=%4.1f ms',carfre/1000,modfre,halflife*1000)); tomwalters@0: tomwalters@0: % plot(sig) tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: return tomwalters@0: end tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: sinus=generatesinus(sig,carfre,amplitude,0); tomwalters@0: tomwalters@0: % calculate envelope and mult both tomwalters@0: envelope=sig; tomwalters@0: tomwalters@0: time_const=halflife/0.69314718; tomwalters@0: tomwalters@0: env_vals=getvalues(envelope); tomwalters@0: sr=getsr(envelope); tomwalters@0: reprate=1/modfre; tomwalters@0: tomwalters@0: tomwalters@0: sig_len=getlength(sig); tomwalters@0: % when regular, all modulataions are at the same time: tomwalters@0: if jitter==0 tomwalters@0: pulse_times=0:reprate:sig_len; tomwalters@0: else tomwalters@0: nr_pulses=ceil(sig_len/reprate); tomwalters@0: modulation_period=1/modfre; tomwalters@0: pulse_times(1)=0; tomwalters@0: for n = 2:nr_pulses tomwalters@0: jittering=(rand-0.5*jitter)*modulation_period; tomwalters@0: pulse_times(n) = pulse_times(n-1) + modulation_period+jittering; tomwalters@0: end tomwalters@0: end tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: % oldval=0; tomwalters@0: % corr=exp(1/sr/time_const); tomwalters@0: % t=0; tomwalters@0: next_pulse=pulse_times(1); tomwalters@0: pulse_counter=1; tomwalters@0: tomwalters@0: % for i=1:getnrpoints(envelope); tomwalters@0: % oldval=oldval/corr; tomwalters@0: % t=t+1/sr; tomwalters@0: % if t>next_pulse tomwalters@0: % oldval=1; tomwalters@0: % pulse_counter=pulse_counter+1; tomwalters@0: % if length(pulse_times)>=pulse_counter tomwalters@0: % next_pulse=pulse_times(pulse_counter); tomwalters@0: % else tomwalters@0: % next_pulse=inf; tomwalters@0: % end tomwalters@0: % end tomwalters@0: % env_vals(i)= oldval; tomwalters@0: % end tomwalters@0: tomwalters@0: % onsettimeconstant tomwalters@0: % tomwalters@0: t=[0:1/sr:reprate-1/sr]*500; tomwalters@0: env=power(t,4).*exp(-2*pi*t/2); tomwalters@0: env=env./max(env); tomwalters@0: tomwalters@0: % e=linspace(env(length(env)-l),0,l+1); tomwalters@0: % env(length(env)-l:end)=e; tomwalters@0: tomwalters@0: % e=linspace(1,0,l+1); tomwalters@0: % env(length(env)-l:end)=env(length(env)-l:end).*e.*e; tomwalters@0: tomwalters@0: t=0; tomwalters@0: ct=1; tomwalters@0: for i=1:getnrpoints(envelope); tomwalters@0: t=t+1/sr; tomwalters@0: ct=ct+1; tomwalters@0: if t>next_pulse tomwalters@0: ct=1; tomwalters@0: pulse_counter=pulse_counter+1; tomwalters@0: if length(pulse_times)>=pulse_counter tomwalters@0: next_pulse=pulse_times(pulse_counter); tomwalters@0: else tomwalters@0: next_pulse=inf; tomwalters@0: end tomwalters@0: end tomwalters@0: if ct>length(env) tomwalters@0: env_vals(i)= 0; tomwalters@0: else tomwalters@0: env_vals(i)= env(ct); tomwalters@0: end tomwalters@0: end tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: % tomwalters@0: % onsettime=0; tomwalters@0: % for i=1:getnrpoints(envelope); tomwalters@0: % t=t+1/sr; tomwalters@0: % tomwalters@0: % % env_vals(i)= exp(-(time)/time_const); tomwalters@0: % env_vals(i)= power(t,onsettime)*exp(-(t)/time_const); tomwalters@0: % t=mod(t,reprate); tomwalters@0: % tomwalters@0: % end tomwalters@0: envelope=setvalues(envelope,env_vals); tomwalters@0: envelope=envelope/max(envelope)*amplitude; tomwalters@0: envelope=setstarttime(envelope,0); tomwalters@0: tomwalters@0: % set the envelope and the amplitude tomwalters@0: sig=sinus*envelope; tomwalters@0: tomwalters@0: sig=setname(sig,sprintf('damped sinusoid %4.2f kHz, Modulation=%4.1f Hz, halflife=%4.1f ms',carfre/1000,modfre,halflife*1000)); tomwalters@0: tomwalters@0: % plot(sig) tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: tomwalters@0: % tomwalters@0: % tomwalters@0: % tomwalters@0: % % from here on old code, might not work! tomwalters@0: % tomwalters@0: % sinus=generatesinus(sig,carfre,amplitude,0); tomwalters@0: % tomwalters@0: % % calculate envelope and mult both tomwalters@0: % envelope=sig; tomwalters@0: % tomwalters@0: % time_const=halflife/0.69314718; tomwalters@0: % tomwalters@0: % env_vals=getvalues(envelope); tomwalters@0: % sr=getsr(envelope); tomwalters@0: % reprate=1/modfre; tomwalters@0: % tomwalters@0: % tomwalters@0: % sig_len=getlength(sig); tomwalters@0: % % when regular, all modulataions are at the same time: tomwalters@0: % if jitter==0 tomwalters@0: % pulse_times=0:reprate:sig_len; tomwalters@0: % else tomwalters@0: % nr_pulses=ceil(sig_len/reprate); tomwalters@0: % modulation_period=1/modfre; tomwalters@0: % pulse_times(1)=0; tomwalters@0: % for n = 2:nr_pulses tomwalters@0: % jittering=(rand-0.5*jitter)*modulation_period; tomwalters@0: % pulse_times(n) = pulse_times(n-1) + modulation_period+jittering; tomwalters@0: % end tomwalters@0: % end tomwalters@0: % tomwalters@0: % tomwalters@0: % tomwalters@0: % % oldval=0; tomwalters@0: % % corr=exp(1/sr/time_const); tomwalters@0: % % t=0; tomwalters@0: % next_pulse=pulse_times(1); tomwalters@0: % pulse_counter=1; tomwalters@0: % tomwalters@0: % % for i=1:getnrpoints(envelope); tomwalters@0: % % oldval=oldval/corr; tomwalters@0: % % t=t+1/sr; tomwalters@0: % % if t>next_pulse tomwalters@0: % % oldval=1; tomwalters@0: % % pulse_counter=pulse_counter+1; tomwalters@0: % % if length(pulse_times)>=pulse_counter tomwalters@0: % % next_pulse=pulse_times(pulse_counter); tomwalters@0: % % else tomwalters@0: % % next_pulse=inf; tomwalters@0: % % end tomwalters@0: % % end tomwalters@0: % % env_vals(i)= oldval; tomwalters@0: % % end tomwalters@0: % tomwalters@0: % % onsettimeconstant tomwalters@0: % % tomwalters@0: % t=[0:1/sr:reprate-1/sr]*500; tomwalters@0: % env=power(t,4).*exp(-2*pi*t/2); tomwalters@0: % env=env./max(env); tomwalters@0: % tomwalters@0: % % e=linspace(env(length(env)-l),0,l+1); tomwalters@0: % % env(length(env)-l:end)=e; tomwalters@0: % tomwalters@0: % % e=linspace(1,0,l+1); tomwalters@0: % % env(length(env)-l:end)=env(length(env)-l:end).*e.*e; tomwalters@0: % tomwalters@0: % t=0; tomwalters@0: % ct=1; tomwalters@0: % for i=1:getnrpoints(envelope); tomwalters@0: % t=t+1/sr; tomwalters@0: % ct=ct+1; tomwalters@0: % if t>next_pulse tomwalters@0: % ct=1; tomwalters@0: % pulse_counter=pulse_counter+1; tomwalters@0: % if length(pulse_times)>=pulse_counter tomwalters@0: % next_pulse=pulse_times(pulse_counter); tomwalters@0: % else tomwalters@0: % next_pulse=inf; tomwalters@0: % end tomwalters@0: % end tomwalters@0: % if ct>length(env) tomwalters@0: % env_vals(i)= 0; tomwalters@0: % else tomwalters@0: % env_vals(i)= env(ct); tomwalters@0: % end tomwalters@0: % end tomwalters@0: % tomwalters@0: % tomwalters@0: % tomwalters@0: % % tomwalters@0: % % onsettime=0; tomwalters@0: % % for i=1:getnrpoints(envelope); tomwalters@0: % % t=t+1/sr; tomwalters@0: % % tomwalters@0: % % % env_vals(i)= exp(-(time)/time_const); tomwalters@0: % % env_vals(i)= power(t,onsettime)*exp(-(t)/time_const); tomwalters@0: % % t=mod(t,reprate); tomwalters@0: % % tomwalters@0: % % end tomwalters@0: % envelope=setvalues(envelope,env_vals); tomwalters@0: % envelope=envelope/max(envelope)*amplitude; tomwalters@0: % envelope=setstarttime(envelope,0); tomwalters@0: % tomwalters@0: % % set the envelope and the amplitude tomwalters@0: % sig=sinus*envelope; tomwalters@0: % tomwalters@0: % sig=setname(sig,sprintf('damped sinusoid %4.2f kHz, Modulation=%4.1f Hz, halflife=%4.1f ms',carfre/1000,modfre,halflife*1000)); tomwalters@0: % tomwalters@0: % % plot(sig)