function [signalout,sf]=supersine(flist,alist,plist,quietflag);
% SUPERSINE Superposition of sines
% function [signalout,sf]=supersine(flist,alist,plist,quietflag)
% supersine: version 30.10.09
%
%   Syntax
%       flist, alist, plist: List of frequencies, amplitudes and phases (the latter in degrees)
%       flist, alist and plist should be row vectors of same length
%       Exception: If alist or plist have only one element then the same amplitude
%       (or phase) is used for all elements of flist
%           Use square brackets for an arbitrary list e.g [100 500 600 800]
%           For a regular series no brackets are necessary e.g
%           100:100:700 corresponds to [100 200 300 400 500 600 700]
%       quietflag: Optional. If present and true, do no sound or graphics
%       signalout and sf: Optional. Returns the synthesized signal and its
%           samplerate
%
%   Hints
%       Some approximations (change n to adjust the number terms as desired):
%           square wave: n=21;supersine((1:2:n)*100,1./(1:2:n),0)
%           sawtooth: n=10;supersine((1:n)*100,(1./(1:n)).*repmat([1 -1],[1 n/2]),0)
%       The effect of using random phase can be examined by using the following expressions for the phase list: 
%           rand(1,ceil(n/2))*360 for the square wave
%           rand(1,n)*360 for the sawtooth wave
%
%   Notes
%       Sound can be heard again by clicking in the figure
%       X axis scaling is set to 50ms, which is intended to show a couple
%       of pitch periods at normal speech fundamental frequencies.
%       The time scale can be changed if desired, e.g to show slow
%       fluctuations in the data when the component sines are not in a
%       simple harmonic relationship, by typing
%       xlim([0 2]) etc.


bequiet=0;
if nargin>3 bequiet=quietflag; end;

if ~bequiet figure; end;


samplerate=10000;
nyquistf=samplerate/2;
sinelength_sec=2;
plotlength=0.050;

sinelength_samples=sinelength_sec*samplerate;


vv=find(flist>=nyquistf);
if ~isempty(vv)
	disp(['Frequency list must not contain values >= ' int2str(nyquistf)]);
	disp('These items have been removed from the list');
	flist(vv)=[];
end;

ns=length(flist);
na=length(alist);
np=length(plist);
if na~=ns 
	if na>1
		disp('Amplitude list of more than one element must match length of frequency list');
		disp('Using first term in list for all frequencies');
	end;
	
	alist=repmat(alist(1),[1 ns]);
	na=1;
end;
if np~=ns 
	if np>1
		disp('Phase list of more than one element must match length of frequency list');
		disp('Using first term in list for all frequencies');
	end;

	plist=repmat(plist(1),[1 ns]);
	np=1;
end;





%for figure title ==============
maxchar=60;
fstr=['f: ' int2str(flist)];
astr=['A: ' num2str(alist(1:na))];
pstr=['Phi: ' int2str(plist(1:np))];

if length(fstr)>maxchar fstr=[fstr(1:maxchar) ' ...']; end;
if length(astr)>maxchar astr=[astr(1:maxchar) ' ...']; end;
if length(pstr)>maxchar pstr=[pstr(1:maxchar) ' ...']; end;
%==================


plistr=plist*(pi/180);
    


mysignal=zeros(sinelength_samples,1);

    tvec=((0:(sinelength_samples-1))/samplerate)';

for ii=1:ns
    
    tmp=alist(ii)*sin((2*pi*flist(ii)*tvec)+plistr(ii));
    mysignal=mysignal+tmp;
end;

if ~bequiet
%soundsc(mysignal,samplerate);
sound((mysignal/max(abs(mysignal)))*0.9,samplerate);
userdata=samplerate;
set(gcf,'userdata',userdata);

hl=plot(tvec,mysignal);
set(gca,'xlim',[0 plotlength]);
grid on;
xlabel('Time (s)','fontsize',14,'fontweight','bold');
ylabel('Ampitude','fontsize',14,'fontweight','bold');

htit=title({fstr,astr,pstr});
set(htit,'vertical','base','fontsize',14,'fontweight','bold');

%activate callback, so sounds can be heard by clicking in any figure
set(gcf,'windowbuttondownfcn',{@playcb,hl});

end;
%keyboard;
if nargout 
    signalout=mysignal; 
    sf=samplerate;
end;
function playcb(obj,data,hl);
samplerate=get(obj,'userdata');
mysignal=get(hl,'ydata');
%soundsc(get(hl,'ydata'),samplerate);
sound((mysignal/max(abs(mysignal)))*0.9,samplerate);
