function out=sampling(action,r)
% SAMPLING  demonstrates unlimited sampling of data, which is
% taken from continuous mathematical functions.
% Push the Return button to return arrays of time, data,
% frequency axis and spectrum to the main workspace.
% SAMPLING('myfun','F') replaces the built-in function SQUARE2
% with the user-defined function F, which is of the form y=F(t).

if nargin==2 & isnumeric(r)		% ODEFile functionality
	rp=[0;0;0];
	rp(1)=-r(2)-r(3);			% Rossler x'
	rp(2)=r(1)+0.15*r(2);		% Rossler y'
	rp(3)=0.2+r(3)*(r(1)-10);	% Rossler z'
	out=rp;

elseif nargin==1 & isstr(action)	% callback functionality
	p=get(gcbf,'UserData');
	switch (action)
	case 'popwave'
		p.wavemode=get(gcbo,'Value');
		if (p.wavemode<3)
			set(p.defb,'Enable','on')
		else
			set(p.defb,'Enable','off')
		end
		p=doit(p);
	case 'define'
		prompt={'Amplitude';'Frequency';'Delay'};
		if p.wavemode==1
			titlestr='Single Sine';
			str=inputdlg(prompt,titlestr,1, ...
				{num2str(p.A),num2str(p.f),num2str(p.delay)});
			if size(str) 
				p.A=eval(str{1});
				p.f=eval(str{2});
				p.delay=eval(str{3});
			end
		else	% defaults are the Fourier series for a square wave
			str=inputdlg('Enter number of periods','P',1,{'3'});
			if size(str)
				p.P=str2num(str{1});
				for k=1:p.P
					titlestr=['Period #',num2str(k)];
					fk=2*k-1;
					Ak=num2str(4/fk/pi);
					fk=num2str(fk);
					str=inputdlg(prompt,titlestr,1,{Ak;fk;'0'});
					if (~size(str)) break, end
					p.Ap(k)=eval(str{1});
					p.fp(k)=eval(str{2});
					p.delayp(k)=eval(str{3});
				end
			end
		end
		p=doit(p);
	case 'fs'
		p.fs=max([eps^0.7,str2num(get(gcbo,'String'))]);
		set(gcbo,'String',num2str(p.fs))
		p=doit(p);
	case 'N'
		p.N=max([1,round(str2num(get(gcbo,'String')))]);
		set(gcbo,'String',num2str(p.N))
		p=doit(p);
	case 'popspect'
		p.spectmode=get(gcbo,'Value');
		p=doit(p);
	case 'rbutton'
		assignin('base','tout',p.t)
		assignin('base','yout',p.y)
		assignin('base','freq',p.freq)
		switch (p.spectmode)
		case 1
			assignin('base','AY',2*abs(p.Y)/p.N)
			disp(['Returned data and Amplitude Spectrum to ', ...
				'[tout,yout] and [freq,AY]'])
		case 2
			assignin('base','PY',2*p.Y.*conj(p.Y)/p.N^2)
			disp(['Returned data and Power Spectrum to ', ...
				'[tout,yout] and [freq,PY]'])
		case 3
			assignin('base','CY',2*p.Y/p.N)
			disp(['Returned data and Complex Spectrum to ', ...
				'[tout,yout] and [freq,CY]'])
		end		
	case 'p'
		out=p;
	case 'help'
    	helpwin('aboutsampling')
	otherwise
		disp(['callback ',action,' unknown'])
	end
	set(gcbf,'UserData',p)

elseif nargin==0 | (nargin==2 & isstr(r))	% start-ups
	myfun=0;	if (nargin==2 & isstr(r)) myfun=1;	end
	p=struct('wavemode',1, ...
		'spectmode',1, ...
		'A',10,'f',1,'delay',0, ...
		'P',3,'Ap',4./[1 3 5]/pi, ...
		'fp',[1 3 5],'delayp',[0 0 0], ...
		'fcnname','square2', ...
		'fs',6, ...
		'N',18, ...
		't',zeros(18,1), ...
		'y',zeros(18,1), ...
		'freq',zeros(9,1), ...
		'Y',zeros(9,1), ...
		'datax',0,'spectax',0,'defb',0);
	if myfun 		% case of custom startup
		p.fcnname=r;
		p.wavemode=4;
	end
    a=figure('Units','normalized','Position',[0.2 0.1 0.75 0.8], ...
      'Numbertitle','off','Name','Sampling','Tag','SamplingBR549');
	set(0,'ShowHiddenHandles','on')
	wh=findobj('Tag','SamplingBR549');
	s=length(wh);
	if (s>1)
		set(a,'Name',['Sampling #',num2str(s)]);
	end
	set(0,'ShowHiddenHandles','off')
	currdir=pwd;
	homedir=pathto('mainswitch.m');
	cd(homedir)
	d=dir('fontpref.mat');
	if (length(d))
		load fontpref
		set(a,'DefaultAxesFontName',axfont.FontName, ...
			'DefaultAxesFontUnits',axfont.FontUnits, ...
			'DefaultAxesFontSize',axfont.FontSize, ...
			'DefaultAxesFontWeight',axfont.FontWeight, ...
			'DefaultAxesFontAngle',axfont.FontAngle, ...
			'DefaultUicontrolFontName',uifont.FontName, ...
			'DefaultUicontrolFontUnits',uifont.FontUnits, ...
			'DefaultUicontrolFontSize',uifont.FontSize, ...
			'DefaultUicontrolFontWeight',uifont.FontWeight, ...
			'DefaultUicontrolFontAngle',uifont.FontAngle)
	end
	cd(currdir)
	color=get(a,'Color');
	p.datax=axes('Parent',a,'Units','normalized','XGrid','on', ...
		'YGrid','on','Position',[0.33 0.57 0.6 0.38]);
	p.spectax=axes('Parent',a,'Units','normalized','XGrid','on', ...
		'YGrid','on','Position',[0.33 0.07 0.6 0.38]);
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.92 0.15 0.05],'String','Signal');
	c=uicontrol('Parent',a,'Style','popupmenu','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.85 0.195 0.05], ...
	   'String',['Sine|Multi-Periodic|Rossler x|',p.fcnname], ...
	   'Callback','sampling popwave','Value',p.wavemode);
	p.defb=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
	   'Position',[0.015 0.770 0.15 0.048],'String','Define', ...
	   'callback','sampling define');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.690 0.15 0.048],'String','Sampling');
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.620 0.15 0.048],'String',num2str(p.fs), ...
	   'callback','sampling fs');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.620 0.045 0.048],'String','fs');
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.540 0.15 0.048],'String',num2str(p.N), ...
	   'callback','sampling N');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.540 0.03 0.048],'String','N');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.3 0.15 0.05],'String','Spectrum');
	c=uicontrol('Parent',a,'Style','popupmenu','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.24 0.195 0.05], ...
	   'String','Amplitude|Power|Complex', ...
	   'Callback','sampling popspect','Value',p.spectmode);
	c=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
	   'Position',[0.015 0.015 0.15 0.05],'String','Return', ...
	   'Callback',['sampling rbutton']);
	b = uimenu('Parent',a,'Label','Sampling');
	c1 = uimenu('Parent',b,'Label','Help','Enable','on','Tag','helpmenu', ...
		'Callback','sampling help');
	c2 = uimenu('Parent',b,'Label','Visible','Enable','on','Tag','onmenu', ...
		'Callback',['set(gcf,''Handlevisibility'',''on'');', ...
			'set(get(gcbo,''UserData''),''Enable'',''on'');', ...
			'set(gcbo,''Enable'',''off'')']);
	c3 = uimenu('Parent',b,'Label','Callback','Enable','off', ...
		'Callback',['set(gcf,''Handlevisibility'',''callback'');', ...
			'set(get(gcbo,''UserData''),''Enable'',''on'');', ...
			'set(gcbo,''Enable'',''off'')'],'UserData',c2);
	set(c2,'UserData',c3)
% 	c=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
% 	   'Position',[0.18 0.07 0.05 0.06],'String','p', ...
% 	   'Callback','p=sampling(''p'')');

	if myfun
		p.wavemode=4;
		set(p.defb,'Enable','off')
	end
	p=doit(p);
	set(a,'HandleVisibility','callback')
	set(a,'UserData',p)

else
	error('Argument Error!')
end

function p=doit(p)	% update displays
p.t=[0:p.N-1]'/p.fs;		% sample time basis
switch (p.wavemode)
case 1		% simple sine
	p.y=p.A*sin(2*pi*p.f*(p.t-p.delay));
	fac=min([floor(1024/p.N),max([2,ceil(20*p.f/p.fs)])]);
	t2=[0:p.N*fac-1]'/(p.fs*fac);	% smooth curve time basis
	y2=p.A*sin(2*pi*p.f*(t2-p.delay));
case 2		% multiperiodic  function
	fac=min([floor(1024/p.N),max([2,ceil(20*max(p.fp)/p.fs)])]);
	t2=[0:p.N*fac-1]'/(p.fs*fac);	% smooth curve time basis
	p.y=zeros(p.N,1);
	y2=zeros(p.N*fac,1);
	for n=1:p.P
		p.y=p.y + p.Ap(n)*sin(2*pi*p.fp(n)*(p.t-p.delayp(n)));
		y2=y2 + p.Ap(n)*sin(2*pi*p.fp(n)*(t2-p.delayp(n)));
	end
case 3		% Rossler x component
	fac=min([floor(1024/p.N),max([2,ceil(20/(2*pi*p.fs))])]);
	t2=[0:p.N*fac-1]'/(p.fs*fac);	% smooth curve time basis
	t2=t2+eps^0.75;		% tiny offset so ode45 won't complain
	tg=[p.t;t2];		% ode45 requires sorted time basis
	[ts,in]=sort(tg);	% set up sort reversal index jn
	[is,jn]=sort(in);
	[td,r]=ode45('Sampling',ts,[0;-10;0]);
	x=r(:,1);
	p.y=x(jn(1:p.N));
	y2=x(jn(p.N+1:p.N*(fac+1)));
case 4		% sample a generic function
	fac=min([floor(1024/p.N),20]);	% typically 20 points/sample
	t2=[0:p.N*fac-1]'/(p.fs*fac);	% smooth curve time basis
	t2=t2+eps^0.75;			% tiny offset in case user function
	tg=[p.t;t2];			% behaves similar to ode45 in disliking
	[ts,in]=sort(tg);		% duplicates in the time basis
	[is,jn]=sort(in);
	y=feval(p.fcnname,ts);
	p.y=y(jn(1:p.N));
	y2=y(jn(p.N+1:p.N*(fac+1)));
end
axes(p.datax)
cla
plot(t2,y2,'r'), hold on
stem(p.t,p.y,'b')
grid on
xlabel('time'),ylabel('y'),title('Function and Samples')
half=floor(p.N/2);
p.freq=[0:half]'/p.N*p.fs;
Y=fft(p.y);
p.Y=Y(1:half+1);
axes(p.spectax)
cla
switch (p.spectmode)
case 1
	stem(p.freq,2*abs(p.Y)/p.N);
	hold on
	switch (p.wavemode)
	case 1
		if p.fs<=p.f*2
			plot(p.f*[1 1],[0 p.A],'r','LineWidth',4)
		end
	case 2
		for k=1:p.P
			if p.fs<=p.fp(k)*2
				plot(p.fp(k)*[1 1],[0 p.Ap(k)],'r','LineWidth',4)
			end
		end		
	case 3
		if p.fs<=1/pi
			plot([1 1]/2/pi,[0 10],'r','LineWidth',4)
		end
	end
	xlabel('Frequency'),ylabel('A')
	title('Amplitude Spectrum')
case 2
	stem(p.freq,2*p.Y.*conj(p.Y)/p.N^2)
	hold on
	switch (p.wavemode)
	case 1
		if p.fs<=p.f*2
			plot(p.f*[1 1],[0 p.A^2/2],'r','LineWidth',4)
		end
	case 2
		for k=1:p.P
			if p.fs<=p.fp(k)*2
				plot(p.fp(k)*[1 1],[0 p.Ap(k)^2/2],'r','LineWidth',4)
			end
		end		
	case 3
		if p.fs<=1/pi
			plot([1 1]/2/pi,[0 50],'r','LineWidth',4)
		end
	end
	xlabel('Frequency'),ylabel('P')
	title('Power Spectrum')
case 3
	plot([1 1]*p.freq(1),[0 real(Y(1))/p.N],'r'), hold on
	plot([1 1]*p.freq(1),[0 imag(Y(1))/p.N],'b')
	stem(p.freq,2*real(p.Y)/p.N,'r')
	stem(p.freq,2*imag(p.Y)/p.N,'b'), grid on
	legend('real','imag',0)
	xlabel('Frequency'),ylabel('Y')
	title('Complex Spectrum')
end
grid on
set(p.spectax,'XLim',[0 p.fs/2])

function s = square2(t)
tmp = mod(t+eps^0.75,2*pi);
nodd = (tmp < pi);
s = 2*nodd-1;

function dest=pathto(fname)
r=which(fname);
[tok1,r]=strtok(r,filesep);
dest=tok1;
[tok2,r]=strtok(r,filesep);
while ~isempty(r)
	tok1=tok2;
	dest=[dest,filesep,tok1];
	[tok2,r]=strtok(r,filesep);
end

% by Tim Conover 1/24/2000
