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