function out=funspect(action)
% FUNSPECT  periodic function spectrum demonstration
% There are no command line arguments, except 'help.'
% Push the Return button to return arrays of time, data,
% frequency axis and spectrum to the main workspace.

if nargin==0		% start-up
	p=struct('wavemode',1, ...
		'spectmode',1, ...
		'A',10, ...
		'f',1, ...
		'delay',0.15, ...
		'S',8, ...
		'M',3, ...
		'fs',6, ...
		'N',18, ...
		't',zeros(18,1), ...
		'y',zeros(18,1), ...
		'freq',zeros(9,1), ...
		'Y',zeros(9,1), ...
		'datax',0,'spectax',0,'trivia',0);
    a=figure('Units','normalized','Position',[0.2 0.1 0.75 0.8], ...
      'Numbertitle','off','Name','FunSpect','Tag','FunSpectBR549');
	set(0,'ShowHiddenHandles','on')
	wh=findobj('Tag','FunSpectBR549');
	s=length(wh);
	if (s>1)
		set(a,'Name',['FunSpect #',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','Wave Form');
	c=uicontrol('Parent',a,'Style','popupmenu','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.85 0.195 0.05], ...
	   'String','Sine|Square|Triangle|Sawtooth', ...
	   'Callback','funspect popwave','Value',p.wavemode);
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.770 0.15 0.048],'String',num2str(p.A), ...
	   'callback','funspect A');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.770 0.03 0.048],'String','A');
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.690 0.15 0.048],'String',num2str(p.f), ...
	   'callback','funspect f');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.690 0.03 0.048],'String','f');
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.620 0.15 0.048],'String',num2str(p.delay), ...
	   'callback','funspect delay');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.620 0.045 0.048],'String','delay');
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.540 0.15 0.048],'String',num2str(p.S), ...
	   'callback','funspect S');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.540 0.03 0.048],'String','S');
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.015 0.460 0.15 0.048],'String',num2str(p.M), ...
	   'callback','funspect M');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.460 0.03 0.048],'String','M');
	p.trivia=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.3 0.22 0.13], ...
	   'String',{'one';'two';'three'});
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.2 0.15 0.05],'String','Spectrum');
	c=uicontrol('Parent',a,'Style','popupmenu','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.015 0.14 0.195 0.05], ...
	   'String','Amplitude|Power|Complex', ...
	   'Callback','funspect popspect','Value',p.spectmode);
	c=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
	   'Position',[0.015 0.015 0.15 0.05],'String','Return', ...
	   'Callback',['funspect rbutton']);
	b = uimenu('Parent',a,'Label','FunSpect');
	c1 = uimenu('Parent',b,'Label','Help','Enable','on','Tag','helpmenu', ...
		'Callback','funspect 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=funspect(''p'')');

	p=doit(p);
	set(a,'HandleVisibility','callback')
	set(a,'UserData',p)

elseif  isstr(action)	% callback functionality
	p=get(gcbf,'UserData');
	switch (action)
	case 'popwave'
		p.wavemode=get(gcbo,'Value');
		p=doit(p);
	case 'popspect'
		p.spectmode=get(gcbo,'Value');
		p=doit(p);
	case {'A','f','delay'}
		eval(['p.',action,'=str2num(get(gcbo,''String''));'])
		p=doit(p);
	case 'S'
		p.S=max([2,round(str2num(get(gcbo,'String')))]);
		set(gcbo,'String',num2str(p.S))
		p=doit(p);		
	case 'M'
		p.M=max([1,round(str2num(get(gcbo,'String')))]);
		set(gcbo,'String',num2str(p.M))
		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('aboutfunspect')
	otherwise
		disp(['callback ',action,' unknown'])
	end
	set(gcbf,'UserData',p)
else
	error('Argument Error!')
end

function p=doit(p)	% update displays
p.fs=p.S*abs(p.f);
p.N=p.S*p.M;
p.t=[0:p.N-1]'/p.fs;
ho=max([1,floor((p.S-1)/2)]);	% harmonic order
fac=min([floor(1024/p.N),max([2,ceil(20*abs(p.f)*ho/p.fs)])]);
t2=[0:p.N*fac-1]'/(p.fs*fac);	% smooth curve time basis
p.y=zeros(p.N,1);
half=floor(p.N/2);
p.freq=[0:half]'/p.N*p.fs;
p.Y=zeros(half+1,1);
A0=0;
An=zeros(ho,1);	Bn=An;
switch (p.wavemode)		% calculate Fourier coefficients
case 1	% sinusoid
	Bn(1)=p.A;
	ho=1;
case 2	% square
	for n=1:2:ho
		Bn(n)=4*p.A/(n*pi);
	end
case 3	% triangle
	for n=1:2:ho
		An(n)=-8*p.A/(n*pi)^2;
	end
case 4	% sawtooth
	for n=1:ho
		Bn(n)=-2*p.A/(n*pi);
		if mod(n,2) Bn(n)=-Bn(n); end
	end
end
p.y=ones(p.N,1)*A0;
y2=ones(p.N*fac,1)*A0;
for n=1:ho			% build  partial Fourier sum and sample it
	wn=2*pi*n*p.f;	% at times t
	p.y=p.y+An(n)*cos(wn*(p.t-p.delay))+Bn(n)*sin(wn*(p.t-p.delay));
	y2=y2+An(n)*cos(wn*(t2-p.delay))+Bn(n)*sin(wn*(t2-p.delay));
end
axes(p.datax)
cla
plot(t2,y2,'r'), hold on
stem(p.t,p.y)
grid on
xlabel('time'),ylabel('y'),title('Periodic Function and Samples')
set(p.trivia,'String',{['harmonics ',num2str(ho)]; ...
	['fs = ',num2str(p.fs)];['N = ',num2str(p.N)]})
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);
	xlabel('Frequency'),ylabel('A')
	title('Amplitude Spectrum')
case 2
	stem(p.freq,2*p.Y.*conj(p.Y)/p.N^2)
	xlabel('Frequency'),ylabel('P')
	title('Power Spectrum')
case 3
	plot([1 1]*p.freq(1),[0 real(Y(1))],'r'), hold on
	plot([1 1]*p.freq(1),[0 imag(Y(1))],'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 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 1999
