Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/~rastor/Software/FourierPeriodogram.m Дата изменения: Tue Feb 18 13:20:41 2014 Дата индексирования: Fri Feb 28 01:21:54 2014 Кодировка: Windows-1251

function [EpochMin,Period,Scatter]=FourierPeriodogram(fname,order,N1st,Pmin,Pmax,delta,epoch);
%
% -------------------------------------------------------------------------
% Program calculates Fourier power spectrum of the time series
%
% Program call:
% [EpochMin,Period,Scatter]=FourierPeriodogram(fname,order,N1st,Pmin,Pmax,delta,epoch)
% or
% FourierPeriodogram(fname,order,N1st,Pmin,Pmax,delta,epoch)
%
% Input:
% fname file name with the time series
% 1st column - time
% 2nd column - data
% order Fourier expansion order
% N1st the number of garmonics (from the 1st one)
% used to estimate the power spectrum by summing
% squared Fourier coefficients
% Pmin,Pmax minimal and maximal trial period values
% delta phase shift between consecutive trial phases
% epoch initial epoch (any)
%
% Output (screen and variables):
% EpochMin epoch of the minimum
% Period period
% Scatter scatter of the phase curve
%
% Programmer: A.S.Rastorguev (Sternberg Astronomical Institute, Lomonosov
% Moscow State University, 7(495)9391616)
% -------------------------------------------------------------------------
close all
data=load(fname);
JD(:,1)=data(:,1);
mag(:,1)=data(:,2);
NN=order;
[nr,nc]=size(data);

% Начальная частота:
omega=1/Pmax;
% Начало счета циклов:
M=0;
while omega <= 1/Pmin,
M=M+1;
om(M,1)=omega;
for i=1:nr
% Вычисление фаз для текущего пробного периода
t(i,1)=omega*(JD(i)-epoch);
t(i,1)=t(i,1)-floor(t(i,1));
end
%% Вычисление коэффициентов Фурье-разложения для текущего периода:
%
% y = x(1)+x(2)*sin (2pi*t*1) + ... + x(NN+1)*sin (2pi*t*N) +
% + x(NN+2)*cos (2pi*t*1) + ...+x(2*NN+1)*cos (2pi*t*N)
%
clear A
A = zeros(nr,2*NN+1);
for i = 1:nr
A(i,1) = 1;
end
for j = 2:(NN+1)
A(:,j) = sin (2*pi*t.*(j-1));
end
for j = (NN+2):(2*NN+1)
A(:,j) = cos (2*pi*t.*(j-NN-1));
end
% Определение массива коэффициентов:
x=A\mag;
% Сумма квадратов всех коэффициентов Фурье-разложения:
sp_all(M,1)=0;
for j=2:(2*NN+1)
sp_all(M,1)=sp_all(M,1)+x(j)^2;
end
% Сумма квадратов амплитуд первых N1st гармоник
sp1st(M,1)=0;
for j=1:N1st
sp1st(M,1)=sp1st(M,1)+x(j+1)^2+x(j+NN+1)^2;
end
% Восстановленное приближение:
z=A*x;
% Рассеяние данных относительно Фурье-разложения:
ff(M,1)=std(mag-z);
% Увеличиваем частоту на delta для следующего цикла:
omega=omega+delta;
end % while
%
% Вывод периодограммы (спектра мощности)
% Размер дисплея:
screen_size=get(0,'MonitorPosition');
% Размер рисунка:
swidth1=screen_size(3)*5/6;
swidth2=screen_size(3)*2/3;
sheight=screen_size(4)*4/5;
h1=figure('Position',[0 screen_size(4)-sheight-100 swidth1 sheight],'NumberTitle','off','Name','Power spectrum, Scatter and 1/Scatter Graphs');
subplot(311)
plot(om,sp1st,'k.');grid;
ylabel('First garmonics');
subplot(312)
plot(om,ff,'k.');grid;
ylabel('Scatter');
subplot(313)
plot(om,1./ff,'k.');grid;
xlabel('Cycles per day');ylabel('1/Scatter');
%
% Определяем пик периодограммы:
[phase_max,I] = min(ff);
P=1/om(I);
% Минимальное рассеяние:
MinScatter=ff(I);
% Рисование фазовой кривой для наилучшего периода:
omega=om(I);
for i=1:nr
% Вычисление фаз для найденного периода
t(i,1)=omega*(JD(i)-epoch);
t(i,1)=t(i,1)-floor(t(i,1));
end
%% Вычисление коэффициентов Фурье-разложения:
%
% y = x(1)+x(2)*sin (2pi*t*1) + ... + x(NN+1)*sin (2pi*t*N) +
% + x(NN+2)*cos (2pi*t*1) + ...+x(2*NN+1)*cos (2pi*t*N)
%
clear A
A = zeros(nr,2*NN+1);
for i = 1:nr
A(i,1) = 1;
end
for j = 2:(NN+1)
A(:,j) = sin (2*pi*t.*(j-1));
end
for j = (NN+2):(2*NN+1)
A(:,j) = cos (2*pi*t.*(j-NN-1));
end
% Определение массива коэффициентов:
x=A\mag;
% Восстановленное приближение:
z=A*x;
%
% Определение момента минимума по сглаженной кривой (разложению):
[phase_max,J] = min(z);
epoch=JD(J);
% Пересчет разложения с новой начальной эпохой (минимума):
for i=1:nr
% Вычисление фаз для найденного периода
t(i,1)=omega*(JD(i)-epoch);
t(i,1)=t(i,1)-floor(t(i,1));
end
%% Вычисление коэффициентов Фурье-разложения:
%
% y = x(1)+x(2)*sin (2pi*t*1) + ... + x(NN+1)*sin (2pi*t*N) +
% + x(NN+2)*cos (2pi*t*1) + ...+x(2*NN+1)*cos (2pi*t*N)
%
clear A
A = zeros(nr,2*NN+1);
for i = 1:nr
A(i,1) = 1;
end
for j = 2:(NN+1)
A(:,j) = sin (2*pi*t.*(j-1));
end
for j = (NN+2):(2*NN+1)
A(:,j) = cos (2*pi*t.*(j-NN-1));
end
% Определение массива коэффициентов:
x=A\mag;
% Восстановленное приближение:
z=A*x;
% Сортировка фаз в порядке возрастания:
[Y,I] = sort(t,'ascend');
%
h2=figure('Position',[0 screen_size(4)-sheight-100 swidth2 sheight],'NumberTitle','off','Name','Light curve');
plot(Y,z(I),'r',Y,mag(I),'k.');grid
xlabel('Phase');ylabel('Magnitude');
%%
disp(' To Period Scatter ')
disp('--------------------------------------------')
disp(sprintf(' %10.6f %10.6f %10.6f',epoch,P,MinScatter));
%
% Вывод результатов - присвоение переменных:
EpochMin=epoch;
Period=P;
Scatter=MinScatter;
%
end % function