function out=FourCoef(action,y)
% FOURCOEF - Illustration of partial Fourier series sums.  The user selects
% periodic function and Fourier series order in a graphic interface.  
% FOURCOEF plots the partial sums for built-in functions
% FOURCOEF(X,Y) uses the function defined by X and Y data in place of the
% default line segments function.  X and Y must be vectors.

if (nargin==1 & isstr(action))	% callback operation
	F=get(gcbf,'UserData');
	switch (action)
	case 'pop'
	   F.wave=get(gcbo,'Value');
	   if (F.wave > 5)
		   set(F.Tclump,'Visible','off')
	   else
		   set(F.Tclump,'Visible','on')
	   end
	   F=doit(F);
   case 'T'
	   F.T=str2num(get(gcbo,'String'));
		F=doit(F);
	case 'N'
		F.N=str2num(get(gcbo,'String'));
		F=doit(F);
	case 'return'
		assignin('base','A0',F.A0)
		assignin('base','An',F.An)
		assignin('base','Bn',F.Bn)
		evalin('base','A0,disp(''[An Bn] = ''),disp([An Bn])')
	case 'F'
		out=F;
	case 'help'
		helpwin('aboutfourcoef')
	otherwise
	   disp(['callback ',action,' unknown'])
	end
	set(gcbf,'UserData',F)

elseif (nargin==0 | nargin==2)
	if nargin==0					% easy invocation
	   x=[-1 1 exp(1) exp(1) pi exp(2)];  
	   y=2*rand(1,6)-1; 
	   wave=1;
	else						% user data invocation
		x=action(:)	;y=y(:);	% make sure they are columns
		if ~(isnumeric(x) & isnumeric(y))
			error('FourCoef Argument Error!')
		end
		if any(diff(x)<0)
			error('all X must be increasing')
		end
		if (length(x)~=length(y))
			error('X and Y must be of equal length')
		end
		wave=8;
	end
	p=struct('N',10, ...
		'wave',wave, ...
		'A',1, ...
		'T',10, ...
		'x',x, ...
		'y',y, ...
		'A0',0, ...
		'An',0, ...
		'Bn',0, ...
		'waveline',0,'sumline',0,'fstxt',0,'Tclump',[0 0]);
	a=figure('Units','normalized','Position',[0.2 0.2 0.75 0.65], ...
		'NumberTitle','off','Name','FourCoef','Tag','FourCoefBR549');
	set(0,'ShowHiddenHandles','on')
	wh=findobj('Tag','FourCoefBR549');
	s=length(wh);
	if (s>1)
		set(a,'Name',['FourCoef #',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');
	b=axes('Parent',a,'Units','normalized','box','on', ...
		'Position',[0.33 0.09 0.6 0.82]);
	c=line('Parent',b,'XData',[],'YData',[],'Color','k','Tag','waveline', ...
	   'LineWidth',1.5);
	p.waveline=c;
	c=line('Parent',b,'XData',[],'YData',[],'Color','b','Tag','sumline', ...
	   'LineWidth',1.5);
	p.sumline=c;
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.02 0.9 0.2 0.05],'String','Select Function');
	c=uicontrol('Parent',a,'Style','popupmenu','Units','normalized', ...
	   'Position',[0.02 0.79 0.2 0.05],'callback','fourcoef pop', ...
	   'String',['odd square|even triangle|odd parabolas|odd triangle|', ...
	   'odd sawtooth|rectified sine|rectified 3 phase|', ...
	   'line segments'],'Value',wave);
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.02 0.65 0.15 0.08],'String',num2str(p.T), ...
	   'callback','fourcoef T');
	p.Tclump(1)=c;
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.65 0.04 0.08],'String','T (sec)');
	p.Tclump(2)=c;
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.02 0.55 0.2 0.05], ...
	   'String',['fundamental = ',num2str(1/p.T),' Hz']);
	p.fstxt=c;
	c=uicontrol('Parent',a,'Style','edit','Units','normalized', ...
	   'Position',[0.02 0.45 0.15 0.08],'String',num2str(p.N), ...
	   'callback','fourcoef N');
	c=uicontrol('Parent',a,'Style','text','Units','normalized', ...
	   'BackgroundColor',color, ...
	   'Position',[0.18 0.45 0.04 0.08],'String','N');
	c=uicontrol('Parent',a,'Style','pushbutton','Units','normalized', ...
	   'Position',[0.03 0.03 0.18 0.1],'String','Return', ...
	   'Callback','fourcoef return');
	b = uimenu('Parent',a,'Label','FourCoef','Tag','uimenu1');
	c1 = uimenu('Parent',b,'Label','Help','Enable','on','Tag','helpmenu', ...
		'Callback','fourcoef 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.03 0.25 0.06 0.1],'String','F' ...
	%    'Callback','F=fourcoef(''F'')');
	
	p=doit(p);	% start-up
	grid on
	if (wave > 5)
		set(p.Tclump,'Visible','off')
	end
	set(a,'HandleVisibility','callback')
	set(a,'UserData',p)

else
	error('FourCoef Argument Error!')
end


function F=doit(F)
np=F.N*32+1;
T=F.T;
xs=[-0.5:1/(np-1):0.5]'*T;	
A0=0;	An=zeros(F.N,1);	Bn=An;
xlabel('x'), ylabel('y'),title('Partial Sum of Fourier Series')
set(gca,'YLimMode','auto')
switch(F.wave)
case 1			% square wave function
	xl=[-0.5 0 0 0.5]*T;	yl=[-1 -1 1 1];
	for n=1:2:F.N
		Bn(n)=4/n/pi;
	end
case 2			% triangle wave function
	xl=[-0.5 0 0.5]*T;	yl=[1 -1 1];
	for n=1:2:F.N
		An(n)=-8/(pi*n)^2;
	end
case 3			% linked parabolas function
	xl=xs;
	mid=floor(np/2);
	yl=[-16*xl(1:mid).^2/T^2-8*xl(1:mid)/T;
		16*xl(mid+1:np).^2/T^2-8*xl(mid+1:np)/T];
	for n=1:2:F.N
		Bn(n)=-32/(pi*n)^3;
	end
case 4			% odd triangle wave function
	xl=[-0.5 -0.25 0.25 0.5]*T;	yl=[0 -1 1 0];
	for n=1:2:F.N
		Bn(n)=8/(pi*n)^2;
		if mod(n,4)==3
			Bn(n)=-Bn(n);
		end
	end
case 5			% sawtooth wave function
	xl=[-0.5 0.5]*T;	yl=[-1 1];
	for n=1:F.N
		Bn(n)=-2/(pi*n);
		if mod(n,2)
			Bn(n)=-Bn(n);
		end
	end
case 6			% rectified sine wave function
	xs=[-0.5:1/(np-1):0.5]'/60;	T=1/60;
	xl=xs;
	mid=floor(np/2);
	A=120*sqrt(2);
	yl=[-A*sin(2*pi*xl(1:mid)/T); A*sin(2*pi*xl(mid+1:np)/T)];
	A0=2*A/pi;
	for n=2:2:F.N
		An(n)=A/pi*(-2/(n-1)+2/(n+1));
	end
	xlabel('Time (sec)'),ylabel('Volts')
	set(gca,'YLim',[0 180])
case 7			% rectified three-phase sine
	xs=[-0.5:1/(np-1):0.5]'/60;	T=1/60;
	xl=xs;
	A=120*sqrt(2);
	yl=A*max(abs([sin(2*pi*60*xl'-pi/3); sin(2*pi*60*xl'); 
		sin(2*pi*60*xl'+pi/3)]))';
	A0=3*A/pi;
	for n=6:6:F.N
		An(n)=3*A/pi*(1/(1-n)+1/(1+n));
	end
	xlabel('Time (sec)'),ylabel('Volts')	
	set(gca,'YLim',[0 180])
case 8			% piecewise linear function
	x=F.x;	y=F.y;
	xl=x;	yl=y;
	L=length(x);
	T=x(L)-x(1);
	nseg=0;
	a=zeros(L-1,1);	b=a;	xi=a;	xf=a;	% set up line segments
	for k=1:L-1
		if x(k+1)>x(k)
			nseg=nseg+1;
			a(nseg)=(y(k+1)-y(k))/(x(k+1)-x(k));
			b(nseg)=y(k)-a(nseg)*x(k);
			xi(nseg)=x(k);	xf(nseg)=x(k+1);
		end
	end
	for k=1:nseg				% calculate mean
		A0=A0+(a(k)*(xf(k)^2-xi(k)^2)/2+b(k)*(xf(k)-xi(k)))/T;
	end
	for n=1:F.N
		wn=2*pi*n/T;
		A=0;	B=0;
		for k=1:nseg		% calculate Fourier coefficients
			A=A+a(k)/(pi*n)*(xf(k)*sin(wn*xf(k))-xi(k)*sin(wn*xi(k)));
			A=A+a(k)*T/2/(pi*n)^2*(cos(wn*xf(k))-cos(wn*xi(k)));
			A=A+b(k)/(pi*n)*(sin(wn*xf(k))-sin(wn*xi(k)));
			B=B-a(k)/(pi*n)*(xf(k)*cos(wn*xf(k))-xi(k)*cos(wn*xi(k)));
			B=B+a(k)*T/2/(pi*n)^2*(sin(wn*xf(k))-sin(wn*xi(k)));
			B=B-b(k)/(pi*n)*(cos(wn*xf(k))-cos(wn*xi(k)));
		end
		An(n)=A;	Bn(n)=B;
	end
	xs=(x(1):T/(np-1):x(L))';	% set up output series
end
series=ones(np,1)*A0;
for n=1:F.N
	wn=2*pi*n/T;
	series=series+An(n)*cos(wn*xs)+Bn(n)*sin(wn*xs);
end
set(F.waveline,'XData',xl,'YData',yl)
set(F.sumline,'XData',xs,'YData',series)
F.A0=A0;	F.An=An;	F.Bn=Bn;
set(F.fstxt,'String',['fundamental = ',num2str(1/T),' Hz'])

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 9/30/99
