lombscargle.m

这是利用MATLAB做L-S谱分析的程序

		sumc=sumc+(c-s)*(c+s);
	end
	wtau=0.5*atan2(2*sumsh,sumc);
	swtau=sin(wtau);
	cwtau=cos(wtau);
	sums=0;
	sumc=0;
	sumsy=0;
	sumcy=0;
	for jval=1:n
		s=wi(jval);
		c=wr(jval);
		ss=s*cwtau-c*swtau;
		cc=c*cwtau+s*swtau;
		sums=sums+ss^2;
		sumc=sumc+cc^2;
		yy=y(jval)-ave;
		sumsy=sumsy+yy*ss;
		sumcy=sumcy+yy*cc;
		wtemp=wr(jval);
		wr(jval)=(wr(jval)*wpr(jval)-wi(jval)*wpi(jval))+wr(jval);
		wi(jval)=(wi(jval)*wpr(jval)+wtemp*wpi(jval))+wi(jval);
	end
	py(ival)=0.5*(sumcy^2/sumc+sumsy^2/sums)/variance;
	
	%WRITE OUTPUT
	freqs(ival,1)=px(ival);
	freqs(ival,2)=py(ival);
	pnow=pnow+1/(ofac*xdif);
end
close(hbar);

pymax=max(py);
jmax=find(py==pymax);

expy=exp(-pymax);

%effm is an estimate of the number of *independent* frequencies
effm=2*nout/ofac;

%The following computation of prob is valid for small values only:
%    prob=effm*expy.
%The following computation of prob is always valid, but presupposes no data clumping:
%    prob=1-(1-expy)^effm.

if ~isempty(effm) & effm~=0
	prob=1-(1-expy)^effm;
	
	fprintf('\npymax = %f',pymax);
	info{lines}=sprintf('pymax = %f',pymax);
	lines=lines+1;
	fprintf('\nfmax = %f',px(jmax));
	info{lines}=sprintf('fmax = %f',px(jmax));
	lines=lines+1;
	fprintf('\neffm = %f',effm);
	info{lines}=sprintf('effm = %f',effm);
	lines=lines+1;
	fprintf('\nexpy = %f',expy);
	info{lines}=sprintf('expy = %f',expy);
	lines=lines+1;
	fprintf('\nnout = %i',nout);
	info{lines}=sprintf('nout = %i',nout);
	lines=lines+1;
	fprintf('\nalpha = %f',prob);
	info{lines}=sprintf('alpha = %f',prob);
	lines=lines+1;
	
	%STORE RESULTS
	results={pymax px(jmax) prob ofac hifac n cal calamp calfreq win wintype nout fhi};
else
	fprintf('period.m: No frequencies to analyze.');
end %if ~isempty(effm) & effm~=0
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function calibrate
global calamp calfreq n x y

msgstring=sprintf('Calibration defaults: Frequency = 0.5 Hz, Amplitude = 1.0.\nDo you want to use these defaults?');
usedefaults=questdlg(msgstring,'Use Defaults?','Yes','No','Yes');
if strcmp(usedefaults,'Yes')
	calfreq=0.5;
	calamp=1.0;
else
	calfreq = getfreq;
	calamp = getamp;
end
%for i=1:n
%   y(i)=y(i)+calamp*sin(calfreq*2*pi*x(i));
%end
y=y+calamp*sin(calfreq*2*pi*x);
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fhi] = getfhi(command_str)
global fhi
fhi=[];

if nargin < 1
	command_str = 'initialize';
end

if ~strcmp(command_str,'initialize')
	handles = get(gcf,'userdata');
	h_edit = handles(1);
end

scrnsize=get(0,'ScreenSize');
posn=[scrnsize(3)/2-150 scrnsize(4)/2-100 300 200];
frameposn=[10 10 posn(3)-20 posn(4)-20];
stringposn=[30 90 posn(3)-60 25];
textposn=[30 120 posn(3)-60 65];

if strcmp(command_str,'initialize')
	h_fig = figure('position',posn,'resize','off','numbertitle','off');
	h_frm = uicontrol(h_fig,'style','frame','position',frameposn);
	h_edit = uicontrol(h_fig,...
		'callback','getfhi(''Answer'');','style','edit','position',stringposn,'userdata',[]);
	uicontrol(h_fig,'style','text','string','fhi: (Hz)','position',textposn);
	handles = [h_edit];
	set(h_fig,'userdata',handles);
	while ~length(get(h_edit,'userdata'))
		drawnow
	end
	fhi = str2num(get(h_edit,'string'));
	if isempty(fhi)
		close(gcf);
		fhi = getfhi;
	end
	close(gcf);
elseif strcmp(command_str,'Answer')
	set(h_edit,'userdata',1);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hifac] = gethifac(command_str)
global hifac
hifac=[];

if nargin < 1
	command_str = 'initialize';
end

if ~strcmp(command_str,'initialize')
	handles = get(gcf,'userdata');
	h_edit = handles(1);
end

scrnsize=get(0,'ScreenSize');
posn=[scrnsize(3)/2-150 scrnsize(4)/2-100 300 200];
frameposn=[10 10 posn(3)-20 posn(4)-20];
stringposn=[30 90 posn(3)-60 25];
textposn=[30 120 posn(3)-60 65];

if strcmp(command_str,'initialize')
	h_fig = figure('position',posn,'resize','off','numbertitle','off');
	h_frm = uicontrol(h_fig,'style','frame','position',frameposn);
	h_edit = uicontrol(h_fig,...
		'callback','gethifac(''Answer'');','style','edit','position',stringposn,'userdata',[],'string','2');
	uicontrol(h_fig,'style','text','string','hifac:','position',textposn,'fontname','Helvetica','fontsize',8);
	handles = [h_edit];
	set(h_fig,'userdata',handles);
	while ~length(get(h_edit,'userdata'))
		drawnow
	end
	hifac = str2num(get(h_edit,'string'));
	if isempty(hifac)
		close(gcf);
		hifac = gethifac;
	end
	close(gcf);
elseif strcmp(command_str,'Answer')
	set(h_edit,'userdata',1);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ofac] = getofac(command_str)
global ofac
ofac=[];

if nargin < 1
	command_str = 'initialize';
end

if ~strcmp(command_str,'initialize')
	handles = get(gcf,'userdata');
	h_edit = handles(1);
end

scrnsize=get(0,'ScreenSize');
posn=[scrnsize(3)/2-150 scrnsize(4)/2-100 300 200];
frameposn=[10 10 posn(3)-20 posn(4)-20];
stringposn=[30 90 posn(3)-60 25];
textposn=[30 120 posn(3)-60 65];

if strcmp(command_str,'initialize')
	h_fig = figure('position',posn,'resize','off','numbertitle','off','name','ofac:');
	h_frm = uicontrol(h_fig,'style','frame','position',frameposn);
	h_edit = uicontrol(h_fig,...
		'callback','getofac(''Answer'');','style','edit','position',stringposn,'userdata',[],'string','4');
	uicontrol(h_fig,'style','text','string','ofac:','position',textposn,'fontname','Helvetica','fontsize',8);
	handles = [h_edit];
	set(h_fig,'userdata',handles);
	while ~length(get(h_edit,'userdata'))
		drawnow
	end
	ofac = str2num(get(h_edit,'string'));
	if isempty(ofac)
		close(gcf);
		ofac = getofac;
	end
	close(gcf);
elseif strcmp(command_str,'Answer')
	set(h_edit,'userdata',1);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dataid] = getdataid(command_str)
if nargin < 1
	command_str = 'initialize';
end
if ~strcmp(command_str,'initialize')
	handles = get(gcf,'userdata');
	h_edit = handles(1);
end
if strcmp(command_str,'initialize')
	h_fig = figure('units','normalized','position',[0.3689 0.3831 0.2604 0.2315],'numbertitle','off',...
		'menubar','none','name','Get Data ID');
	h_frm = uicontrol(h_fig,'style','frame','units','normalized','position',[0.05 0.05 0.9 0.9]);
	h_edit = uicontrol(h_fig,...
		'callback','getdataid(''Answer'');','style','edit','units', 'normalized', 'position',[0.1 0.4 0.8 0.15],'userdata',[]);
	uicontrol(h_fig,'style','text','string','Enter a data-set identifier:','units','normalized','position',[0.1 0.7 0.8 0.15],...
		'fontname','Helvetica','fontunits', 'normalized', 'fontsize',0.5,'horizontalalignment','center');
	handles = [h_edit];
	set(h_fig,'userdata',handles);
	while ~length(get(h_edit,'userdata'))
		drawnow
	end
	dataid = get(h_edit,'string');
	close(gcf);
elseif strcmp(command_str,'Answer')
	set(h_edit,'userdata',1);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [calfreq] = getfreq(command_str)

if nargin < 1
	command_str = 'initialize';
end

if ~strcmp(command_str,'initialize')
	handles = get(gcf,'userdata');
	h_edit = handles(1);
end

scrnsize=get(0,'ScreenSize');
posn=[scrnsize(3)/2-150 scrnsize(4)/2-100 300 200];
frameposn=[10 10 posn(3)-20 posn(4)-20];
stringposn=[30 90 posn(3)-60 25];
textposn=[30 120 posn(3)-60 65];

if strcmp(command_str,'initialize')
	h_fig = figure('position',posn,'resize','off','numbertitle','off');
	h_frm = uicontrol(h_fig,'style','frame','position',frameposn);
	h_edit = uicontrol(h_fig,...
		'callback','getfreq(''Answer'');','style','edit','position',stringposn,'userdata',[]);
	uicontrol(h_fig,'style','text','string','Calibration Frequency (Hz):','position',textposn,'fontname','Helvetica','fontsize',8);
	handles = [h_edit];
	set(h_fig,'userdata',handles);
	while ~length(get(h_edit,'userdata'))
		drawnow
   end
   calfreq = str2num(get(h_edit,'string'));
	close(gcf);
elseif strcmp(command_str,'Answer')
	set(h_edit,'userdata',1);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [calamp] = getamp(command_str)

if nargin < 1
	command_str = 'initialize';
end

if ~strcmp(command_str,'initialize')
	handles = get(gcf,'userdata');
	h_edit = handles(1);
end

scrnsize=get(0,'ScreenSize');
posn=[scrnsize(3)/2-150 scrnsize(4)/2-100 300 200];
frameposn=[10 10 posn(3)-20 posn(4)-20];
stringposn=[30 90 posn(3)-60 25];
textposn=[30 120 posn(3)-60 65];

if strcmp(command_str,'initialize')
	h_fig = figure('position',posn,'resize','off','numbertitle','off');
	h_frm = uicontrol(h_fig,'style','frame','position',frameposn);
	h_edit = uicontrol(h_fig,...
		'callback','getamp(''Answer'');','style','edit','position',stringposn,'userdata',[]);
	uicontrol(h_fig,'style','text','string','Calibration Amplitude:','position',textposn,'fontname','Helvetica','fontsize',8);
	handles = [h_edit];
	set(h_fig,'userdata',handles);
	while ~length(get(h_edit,'userdata'))
		drawnow
   end
   calamp = str2num(get(h_edit,'string'));
	close(gcf);
elseif strcmp(command_str,'Answer')
	set(h_edit,'userdata',1);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%