view aim-mat/tools/@signal/genharmonics.asv @ 4:537f939baef0 tip

various bug fixes and changed copyright message

% method of class @signal
% function sig=genharmonicstim(sig,varargin)
%   INPUT VALUES:
%       sig: @signal with length and samplerate 
%       varargin must have several parameters:
%       fundamental (default 128) = periodicity
%       min_fre (128) = lowest possible frequency
%       max_fre     (10000) = highest possible frequency
%
% 		envelope of amplitudes
%		either
%       filterprop ([256,256,1024,512]) = fc, df1, bandwidth, df2 (in Hz)
%       default values: fc=3500;
%						df1=256; 
%                       bandwidth=1024;
%                       df2=512;
%		the amplitdes can also be given by cf and two slopes for
%		higher and lower frequencies:
%					eg 'cf',1000,'lowslope',25,'highslope',38
% 			the highest and lowest possible allowed harmonic are given
% 			in either case by giving 'lowestharmonic' and
% 			'highestharmonic' (default value: 1 and inf)
%
%       type =              which type   (default none) 
%                           filtered, 
%                           decreaseoddamplitude, 
%                           decreaseoddphase
%                           shiftallcomponents
%                           mistunedharmonic
%							decrease_amplitude_linear
%       changeby = value, that the odd harmonics shall vary (degree or dB or whatever)
%       phases must be in degrees!
%   RETURN VALUE:
%       sig:  @signal 
%
% examples:	
% create a stimulus with certain filtercharacteristic with random
% phase, and every second harmonic reduced by an amount
% tone(i)=genharmonics(sig,'fundamental',chroma,...
% 		'filterprop',[toneheight,handles.df1,handles.bw,handles.df2],...
% 		'phase','random',...
% 		'type','decreaseoddamplitude',...
% 		'changeby',octheight...
% 	);
% 

% (c) 2003, University of Cambridge, Medical Research Council 
% Stefan Bleeck ([e-mail])
% http://www.mrc-cbu.cam.ac.uk/cnbh/aimmanual
% $Date: 2003/09/19 09:32:22 $
% $Revision: 1.12 $

function sig=genharmonics(sig,varargin)

if mod(nargin,2)==0
    disp('odd number of parameters - please input a full set of parameters and arguments');
    return;
end
str_fundamental=getargument(varargin,'fundamental');
str_type=getargument(varargin,'type');
str_harmnr=getargument(varargin,'harmonicnumber');
str_changeby=getargument(varargin,'changeby');
str_filterprop=getargument(varargin,'filterprop');
str_fc=getargument(varargin,'fc');
str_lowslope=getargument(varargin,'lowslope');
str_highslope=getargument(varargin,'highslope');
str_lowestharmonic=getargument(varargin,'lowestharmonic');
str_highestharmonic=getargument(varargin,'highestharmonic');
str_bw=getargument(varargin,'bw');
str_phase=getargument(varargin,'phase');
str_min_fre=getargument(varargin,'min_fre');
str_max_fre=getargument(varargin,'max_fre');

str_which_harmonics=getargument(varargin,'which harmonics');

% defaultvalues:
if strcmp(str_filterprop,'')
	if strcmp(str_changeby,'')   
		str_filterprop{1}=[256,256,1024,512];
		fc=256;
		df1=256;
		bandwidth=1024;
		df2=512;
	else
		min_fre=str_min_fre{1};
		max_fre=str_max_fre{1};
		df1=1;
		df2=1;
		fc=min_fre;
		bandwidth=max_fre-min_fre;
	end
else
	%eval(sprintf('filterprop=%s;',str_filterprop{1}));
	fc=str_filterprop{1}(1);
	df1=str_filterprop{1}(2);
	bandwidth=str_filterprop{1}(3);
	df2=str_filterprop{1}(4);
end
if strcmp(str_changeby,'')   
else
	if isnumeric(str_changeby{1})
		changeby=str_changeby{1};
	else
		eval(sprintf('changeby=%f;',str_changeby{1}));
	end
end

% different method of defining envelope: center frequency is the
% highest point, highslope and lowslope define the amplitude on both
% sides for each harmonic
if ~strcmp(str_lowslope,'') && ~strcmp(str_highslope,'')
	lowslope=str_lowslope{1};
	highslope=str_highslope{1};
	calculate_amplitude_with_slopes=1;
	
	% test with 
	% plot(powerspectrum(genharmonics(signal(0.1,16000),'fc',2000,'lowslope',30,'highslope',40,'fundamental',250)))
else
	calculate_amplitude_with_slopes=0;
end

if strcmp(str_type,'')    
    str_type{1}='';
    type='';
else
    type=str_type{1};
end

if strcmp(str_phase,'')
	setphase='cosine';
else
    setphase=str_phase{1};
end

if strcmp(str_fundamental,'') 
    str_fundamental{1}='128';
    fundamental=128;
else
	if isnumeric(str_fundamental{1})
		fundamental=str_fundamental{1};
	else
	    eval(sprintf('fundamental=%s;',str_fundamental{1}));
	end
end

if strcmp(str_lowestharmonic,'') 
    lowestharmonic=1;
else
	if isnumeric(str_lowestharmonic{1})
		lowestharmonic=str_lowestharmonic{1};
	else
	    eval(sprintf('lowestharmonic=%s;',lowestharmonic{1}));
	end
end
if strcmp(str_highestharmonic,'') 
    highestharmonic=inf;
else
	if isnumeric(str_highestharmonic{1})
		highestharmonic=str_highestharmonic{1};
	else
	    eval(sprintf('highestharmonic=%s;',highestharmonic{1}));
	end
end

if strcmp(str_fc,'') 
else
	if isnumeric(str_fc{1})
		fc=str_fc{1};
	else
		eval(sprintf('fc=%s;',str_fc{1}));
	end
end

if strcmp(str_bw,'') 
else
	if isnumeric(str_bw{1})
		bandwidth=str_bw{1};
	else
	    eval(sprintf('bandwidth=%s;',str_bw{1}));
	end
end

if strcmp(str_which_harmonics,'') 
    str_which_harmonics{1}='all';
end

samplerate=sig.samplerate;
length=getlength(sig);


% begin!

max_fre=fc+bandwidth+df2;
min_fre=fc-df1;
if (min_fre<0)  %squeese df1 to go from 0 to fc
    df1=df1-abs(min_fre);
%     disp('df1 was reduced')
end

if fundamental > max_fre
    disp('error: genharmonics: fundamental must be smaller then highest frequency');
    return;
end

s=signal(length,samplerate);
if max_fre>1000
	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));
else
	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));
end

if calculate_amplitude_with_slopes
	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));
end

% in case of sloped amplitudes, we dont want a limit on harmonics
if calculate_amplitude_with_slopes
	max_fre=getsr(sig)/2;
	min_fre=0;
end
	
% if limit of harmonics is explicitly given
if ~strcmp(str_highestharmonic,'') 
	max_fre=highestharmonic*fundamental;
end
if ~strcmp(str_lowestharmonic,'') 
	min_fre=lowestharmonic*fundamental;
end
	
fre=fundamental;
count_partials=1;
while fre <= max_fre
    if fre >= min_fre
        temp=signal(length,samplerate);
        amplitude=1;
        phase=0;
        offset=0;
        if strcmp(type,'mistunedharmonic') % in %
            eval(sprintf('nr=%s;',str_harmnr{1}));
            if count_partials==nr
                offset=fundamental*changeby/100;
                amplitude=1;
                phase=0;
            end     
        end
        if strcmp(type,'shiftallcomponents')   % in Hz
            offset=changeby;
            amplitude=1;
            phase=0;
        end
        if strcmp(type,'decreaseoddphase')
            % phase must be given in degree!
            if mod(count_partials,2)==1
                amplitude=1;
                phase=changeby;
            end
        end

		if strcmp(type,'decrease_amplitude_linear')
			% the amount must be given in changeby!
            amplitude=1*power(10,(-changeby*count_partials)/20);
            ampscale=1;
        else
            ampscale=filterbandamp(fre+offset,amplitude,fc,df1,bandwidth,df2);
        end
        amplitude=amplitude*ampscale;


		% stattdessen mit Slopes:
		if calculate_amplitude_with_slopes
			% calculate the distance from cf (in octaves) 
			% and from this the attenuation
			distance=log2(fre/fc);
			if distance >= 0
				atten=distance*highslope;
			else
				atten=-distance*lowslope;
			end
            amplitude=1*power(10,-atten/20);
		end

        if strcmp(type,'decreaseoddamplitude')
            if mod(count_partials,2)==1
                amplitude=amplitude*power(10,changeby/20);
            end
        end
        if strcmp(type,'decreaseevenamplitude')
            if mod(count_partials,2)==0
                amplitude=amplitude*power(10,changeby/20);
            end
        end
        
        
        switch str_which_harmonics{1}
            case 'all'
            case 'only odd'
                if mod(count_partials,2)==0
                    amplitude=0;
                end
            case 'only even'
                if mod(count_partials,2)==1
                    amplitude=0;
                end
        end        
		
        % degree2rad
        switch setphase
            case 'random'
                piphase=rand(1)*2*pi+pi;
            case 'cos'
% 			piphase=phase*(pi/180)+pi/2;
			piphase=phase*(pi/180);
        end
%         disp(sprintf('fre: %3.2f amp:%2.1f',fre,amplitude*100));
%         amplitude
%         fre
        temp=generatesinus(temp,fre+offset,amplitude,piphase);  
        
        % add them up!
        s=s+temp;
    end
    fre=fre+fundamental;
    count_partials=count_partials+1;
end

sig=s;