tomwalters@0: % method of class @signal
tomwalters@0: % function sig=genharmonicstim(sig,varargin)
tomwalters@0: %   INPUT VALUES:
tomwalters@0: %       sig: @signal with length and samplerate 
tomwalters@0: %       varargin must have several parameters:
tomwalters@0: %       fundamental (default 128) = periodicity
tomwalters@0: %       min_fre (128) = lowest possible frequency
tomwalters@0: %       max_fre     (10000) = highest possible frequency
tomwalters@0: %
tomwalters@0: % 		envelope of amplitudes
tomwalters@0: %		either
tomwalters@0: %       filterprop ([256,256,1024,512]) = fc, df1, bandwidth, df2 (in Hz)
tomwalters@0: %       default values: fc=3500;
tomwalters@0: %						df1=256; 
tomwalters@0: %                       bandwidth=1024;
tomwalters@0: %                       df2=512;
tomwalters@0: %		the amplitdes can also be given by cf and two slopes for
tomwalters@0: %		higher and lower frequencies:
tomwalters@0: %					eg 'cf',1000,'lowslope',25,'highslope',38
tomwalters@0: % 			the highest and lowest possible allowed harmonic are given
tomwalters@0: % 			in either case by giving 'lowestharmonic' and
tomwalters@0: % 			'highestharmonic' (default value: 1 and inf)
tomwalters@0: %
tomwalters@0: %       type =              which type   (default none) 
tomwalters@0: %                           filtered, 
tomwalters@0: %                           decreaseoddamplitude, 
tomwalters@0: %                           decreaseoddphase
tomwalters@0: %                           shiftallcomponents
tomwalters@0: %                           mistunedharmonic
tomwalters@0: %							decrease_amplitude_linear
tomwalters@0: %       changeby = value, that the odd harmonics shall vary (degree or dB or whatever)
tomwalters@0: %       phases must be in degrees!
tomwalters@0: %   RETURN VALUE:
tomwalters@0: %       sig:  @signal 
tomwalters@0: %
tomwalters@0: % examples:	
tomwalters@0: % create a stimulus with certain filtercharacteristic with random
tomwalters@0: % phase, and every second harmonic reduced by an amount
tomwalters@0: % tone(i)=genharmonics(sig,'fundamental',chroma,...
tomwalters@0: % 		'filterprop',[toneheight,handles.df1,handles.bw,handles.df2],...
tomwalters@0: % 		'phase','random',...
tomwalters@0: % 		'type','decreaseoddamplitude',...
tomwalters@0: % 		'changeby',octheight...
tomwalters@0: % 	);
tomwalters@0: % 
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/09/19 09:32:22 $
tomwalters@0: % $Revision: 1.12 $
tomwalters@0: 
tomwalters@0: function sig=genharmonics(sig,varargin)
tomwalters@0: 
tomwalters@0: if mod(nargin,2)==0
tomwalters@0:     disp('odd number of parameters - please input a full set of parameters and arguments');
tomwalters@0:     return;
tomwalters@0: end
tomwalters@0: str_fundamental=getargument(varargin,'fundamental');
tomwalters@0: str_type=getargument(varargin,'type');
tomwalters@0: str_harmnr=getargument(varargin,'harmonicnumber');
tomwalters@0: str_changeby=getargument(varargin,'changeby');
tomwalters@0: str_filterprop=getargument(varargin,'filterprop');
tomwalters@0: str_fc=getargument(varargin,'fc');
tomwalters@0: str_lowslope=getargument(varargin,'lowslope');
tomwalters@0: str_highslope=getargument(varargin,'highslope');
tomwalters@0: str_lowestharmonic=getargument(varargin,'lowestharmonic');
tomwalters@0: str_highestharmonic=getargument(varargin,'highestharmonic');
tomwalters@0: str_bw=getargument(varargin,'bw');
tomwalters@0: str_phase=getargument(varargin,'phase');
tomwalters@0: str_min_fre=getargument(varargin,'min_fre');
tomwalters@0: str_max_fre=getargument(varargin,'max_fre');
tomwalters@0: 
tomwalters@0: str_which_harmonics=getargument(varargin,'which harmonics');
tomwalters@0: 
tomwalters@0: % defaultvalues:
tomwalters@0: if strcmp(str_filterprop,'')
tomwalters@0: 	if strcmp(str_changeby,'')   
tomwalters@0: 		str_filterprop{1}=[256,256,1024,512];
tomwalters@0: 		fc=256;
tomwalters@0: 		df1=256;
tomwalters@0: 		bandwidth=1024;
tomwalters@0: 		df2=512;
tomwalters@0: 	else
tomwalters@0: 		min_fre=str_min_fre{1};
tomwalters@0: 		max_fre=str_max_fre{1};
tomwalters@0: 		df1=1;
tomwalters@0: 		df2=1;
tomwalters@0: 		fc=min_fre;
tomwalters@0: 		bandwidth=max_fre-min_fre;
tomwalters@0: 	end
tomwalters@0: else
tomwalters@0: 	%eval(sprintf('filterprop=%s;',str_filterprop{1}));
tomwalters@0: 	fc=str_filterprop{1}(1);
tomwalters@0: 	df1=str_filterprop{1}(2);
tomwalters@0: 	bandwidth=str_filterprop{1}(3);
tomwalters@0: 	df2=str_filterprop{1}(4);
tomwalters@0: end
tomwalters@0: if strcmp(str_changeby,'')   
tomwalters@0: else
tomwalters@0: 	if isnumeric(str_changeby{1})
tomwalters@0: 		changeby=str_changeby{1};
tomwalters@0: 	else
tomwalters@0: 		eval(sprintf('changeby=%f;',str_changeby{1}));
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % different method of defining envelope: center frequency is the
tomwalters@0: % highest point, highslope and lowslope define the amplitude on both
tomwalters@0: % sides for each harmonic
tomwalters@0: if ~strcmp(str_lowslope,'') && ~strcmp(str_highslope,'')
tomwalters@0: 	lowslope=str_lowslope{1};
tomwalters@0: 	highslope=str_highslope{1};
tomwalters@0: 	calculate_amplitude_with_slopes=1;
tomwalters@0: 	
tomwalters@0: 	% test with 
tomwalters@0: 	% plot(powerspectrum(genharmonics(signal(0.1,16000),'fc',2000,'lowslope',30,'highslope',40,'fundamental',250)))
tomwalters@0: else
tomwalters@0: 	calculate_amplitude_with_slopes=0;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if strcmp(str_type,'')    
tomwalters@0:     str_type{1}='';
tomwalters@0:     type='';
tomwalters@0: else
tomwalters@0:     type=str_type{1};
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if strcmp(str_phase,'')
tomwalters@0: 	setphase='cosine';
tomwalters@0: else
tomwalters@0:     setphase=str_phase{1};
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if strcmp(str_fundamental,'') 
tomwalters@0:     str_fundamental{1}='128';
tomwalters@0:     fundamental=128;
tomwalters@0: else
tomwalters@0: 	if isnumeric(str_fundamental{1})
tomwalters@0: 		fundamental=str_fundamental{1};
tomwalters@0: 	else
tomwalters@0: 	    eval(sprintf('fundamental=%s;',str_fundamental{1}));
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if strcmp(str_lowestharmonic,'') 
tomwalters@0:     lowestharmonic=1;
tomwalters@0: else
tomwalters@0: 	if isnumeric(str_lowestharmonic{1})
tomwalters@0: 		lowestharmonic=str_lowestharmonic{1};
tomwalters@0: 	else
tomwalters@0: 	    eval(sprintf('lowestharmonic=%s;',lowestharmonic{1}));
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: if strcmp(str_highestharmonic,'') 
tomwalters@0:     highestharmonic=inf;
tomwalters@0: else
tomwalters@0: 	if isnumeric(str_highestharmonic{1})
tomwalters@0: 		highestharmonic=str_highestharmonic{1};
tomwalters@0: 	else
tomwalters@0: 	    eval(sprintf('highestharmonic=%s;',highestharmonic{1}));
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if strcmp(str_fc,'') 
tomwalters@0: else
tomwalters@0: 	if isnumeric(str_fc{1})
tomwalters@0: 		fc=str_fc{1};
tomwalters@0: 	else
tomwalters@0: 		eval(sprintf('fc=%s;',str_fc{1}));
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if strcmp(str_bw,'') 
tomwalters@0: else
tomwalters@0: 	if isnumeric(str_bw{1})
tomwalters@0: 		bandwidth=str_bw{1};
tomwalters@0: 	else
tomwalters@0: 	    eval(sprintf('bandwidth=%s;',str_bw{1}));
tomwalters@0: 	end
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if strcmp(str_which_harmonics,'') 
tomwalters@0:     str_which_harmonics{1}='all';
tomwalters@0: end
tomwalters@0: 
tomwalters@0: samplerate=sig.samplerate;
tomwalters@0: length=getlength(sig);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % begin!
tomwalters@0: 
tomwalters@0: max_fre=fc+bandwidth+df2;
tomwalters@0: min_fre=fc-df1;
tomwalters@0: if (min_fre<0)  %squeese df1 to go from 0 to fc
tomwalters@0:     df1=df1-abs(min_fre);
tomwalters@0: %     disp('df1 was reduced')
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if fundamental > max_fre
tomwalters@0:     disp('error: genharmonics: fundamental must be smaller then highest frequency');
tomwalters@0:     return;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: s=signal(length,samplerate);
tomwalters@0: if max_fre>1000
tomwalters@0: 	s=setname(s,sprintf('Harmonic Signal - f0=%4.1f Hz, type: %s (from %2.2f kHz to %2.2f kHz)',fundamental,type,min_fre/1000,max_fre/1000));
tomwalters@0: else
tomwalters@0: 	s=setname(s,sprintf('Harmonic Signal - f0=%4.1f Hz, type: %s  (from %3.0 Hz to %3.0f Hz)',fundamental,type,min_fre,max_fre));
tomwalters@0: end
tomwalters@0: 
tomwalters@0: if calculate_amplitude_with_slopes
tomwalters@0: 	s=setname(s,sprintf('Harmonic Signal - modfre=%4.1f Hz, type: %s  (cf: %2.2f kHz, low slope: %3.0f dB/oct, high slope %3.0f dB/oct)',fundamental,type,fc/1000,lowslope,highslope));
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % in case of sloped amplitudes, we dont want a limit on harmonics
tomwalters@0: if calculate_amplitude_with_slopes
tomwalters@0: 	max_fre=getsr(sig)/2;
tomwalters@0: 	min_fre=0;
tomwalters@0: end
tomwalters@0: 	
tomwalters@0: % if limit of harmonics is explicitly given
tomwalters@0: if ~strcmp(str_highestharmonic,'') 
tomwalters@0: 	max_fre=highestharmonic*fundamental;
tomwalters@0: end
tomwalters@0: if ~strcmp(str_lowestharmonic,'') 
tomwalters@0: 	min_fre=lowestharmonic*fundamental;
tomwalters@0: end
tomwalters@0: 	
tomwalters@0: fre=fundamental;
tomwalters@0: count_partials=1;
tomwalters@0: while fre <= max_fre
tomwalters@0:     if fre >= min_fre
tomwalters@0:         temp=signal(length,samplerate);
tomwalters@0:         amplitude=1;
tomwalters@0:         phase=0;
tomwalters@0:         offset=0;
tomwalters@0:         if strcmp(type,'mistunedharmonic') % in %
tomwalters@0:             eval(sprintf('nr=%s;',str_harmnr{1}));
tomwalters@0:             if count_partials==nr
tomwalters@0:                 offset=fundamental*changeby/100;
tomwalters@0:                 amplitude=1;
tomwalters@0:                 phase=0;
tomwalters@0:             end     
tomwalters@0:         end
tomwalters@0:         if strcmp(type,'shiftallcomponents')   % in Hz
tomwalters@0:             offset=changeby;
tomwalters@0:             amplitude=1;
tomwalters@0:             phase=0;
tomwalters@0:         end
tomwalters@0:         if strcmp(type,'decreaseoddphase')
tomwalters@0:             % phase must be given in degree!
tomwalters@0:             if mod(count_partials,2)==1
tomwalters@0:                 amplitude=1;
tomwalters@0:                 phase=changeby;
tomwalters@0:             end
tomwalters@0:         end
tomwalters@0: 
tomwalters@0: 		if strcmp(type,'decrease_amplitude_linear')
tomwalters@0: 			% the amount must be given in changeby!
tomwalters@0:             amplitude=1*power(10,(-changeby*count_partials)/20);
tomwalters@0:             ampscale=1;
tomwalters@0:         else
tomwalters@0:             ampscale=filterbandamp(fre+offset,amplitude,fc,df1,bandwidth,df2);
tomwalters@0:         end
tomwalters@0:         amplitude=amplitude*ampscale;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 		% stattdessen mit Slopes:
tomwalters@0: 		if calculate_amplitude_with_slopes
tomwalters@0: 			% calculate the distance from cf (in octaves) 
tomwalters@0: 			% and from this the attenuation
tomwalters@0: 			distance=log2(fre/fc);
tomwalters@0: 			if distance >= 0
tomwalters@0: 				atten=distance*highslope;
tomwalters@0: 			else
tomwalters@0: 				atten=-distance*lowslope;
tomwalters@0: 			end
tomwalters@0:             amplitude=1*power(10,-atten/20);
tomwalters@0: 		end
tomwalters@0: 
tomwalters@0:         if strcmp(type,'decreaseoddamplitude')
tomwalters@0:             if mod(count_partials,2)==1
tomwalters@0:                 amplitude=amplitude*power(10,changeby/20);
tomwalters@0:             end
tomwalters@0:         end
tomwalters@0:         if strcmp(type,'decreaseevenamplitude')
tomwalters@0:             if mod(count_partials,2)==0
tomwalters@0:                 amplitude=amplitude*power(10,changeby/20);
tomwalters@0:             end
tomwalters@0:         end
tomwalters@0:         
tomwalters@0:         
tomwalters@0:         switch str_which_harmonics{1}
tomwalters@0:             case 'all'
tomwalters@0:             case 'only odd'
tomwalters@0:                 if mod(count_partials,2)==0
tomwalters@0:                     amplitude=0;
tomwalters@0:                 end
tomwalters@0:             case 'only even'
tomwalters@0:                 if mod(count_partials,2)==1
tomwalters@0:                     amplitude=0;
tomwalters@0:                 end
tomwalters@0:         end        
tomwalters@0: 		
tomwalters@0:         % degree2rad
tomwalters@0:         switch setphase
tomwalters@0:             case 'random'
tomwalters@0:                 piphase=rand(1)*2*pi+pi;
tomwalters@0:             case 'cos'
tomwalters@0: % 			piphase=phase*(pi/180)+pi/2;
tomwalters@0: 			piphase=phase*(pi/180);
tomwalters@0:         end
tomwalters@0: %         disp(sprintf('fre: %3.2f amp:%2.1f',fre,amplitude*100));
tomwalters@0: %         amplitude
tomwalters@0: %         fre
tomwalters@0:         temp=generatesinus(temp,fre+offset,amplitude,piphase);  
tomwalters@0:         
tomwalters@0:         % add them up!
tomwalters@0:         s=s+temp;
tomwalters@0:     end
tomwalters@0:     fre=fre+fundamental;
tomwalters@0:     count_partials=count_partials+1;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: sig=s;