spectrum.m

功能： 功率谱分析 算法参考文献：黄嘉佑

function spectrum()
%spectrum()
%%load data and some parameters.
inputfile=input('Please input the data file name:','s');
fidin=fopen(inputfile,'r');
helpin=strcat('Can''t open the file---',inputfile);
if fidin==-1
    error(helpin);
end
x=fscanf(fidin,'%f');
fclose(fidin);
N=length(x);
M=input('Please input the M value:');
alpha=input('Please input the confidence interval(in decimal):');
%Calculating auto_connection coefficient.
aver_x=sum(x)./N;
fangcha_x=sum((x-aver_x).^2)./N;
for c1=0:M
    CC(c1+1)=0;
    for c2=1:(N-c1)
        CC(c1+1)=CC(c1+1)+(x(c2)-aver_x).*(x(c2+c1)-aver_x)./fangcha_x;
    end
    CC(c1+1)=CC(c1+1)./(N-c1);
end
%Calculatting the smoothed power spectra.
for c1=0:M
    if c1==0||c1==M
        Bk=0.5;
    else
        Bk=1;
    end
    right2=0;
    for c2=1:(M-1)
        right2=right2+CC(c2+1).*cos(pi.*c1.*c2./M).*(1+cos(pi.*c2./M));
    end
    sk(c1+1)=Bk.*(right2+CC(1))./M;
end
%Noice type examination.
xindu=1-alpha;
stdcor_x=stdcor(xindu,N);
if CC(2)>=stdcor_x
    noicetype=1; %Red noice.
else
    noicetype=2; %White noice.
end
%Calculatting the stdandard spectra.
v=(2.*N-M./2)./M;
free_x2=chi2inv(alpha,v)./v;
aver_sk=sum(sk)./(M+1);
if noicetype==1 %Red noice.
    for c1=0:M
        standard(c1+1)=free_x2.*aver_sk.*(1-CC(2).^2)./(1+CC(2).^2-2.*CC(2).*cos(pi.*c1./M));
    end
end
if noicetype==2 %White noice.
    standard(c1+1)=aver_sk.*free_x2;
end
%Calculate the length of cycle.
T=0:M;
T=2.*N./T;
axis_x=0:M;
%output results.
outputfile=input('Please input the outputfile name:','s');
fidout1=fopen(outputfile,'w');
helpword=strcat('Can''t open the file----',outputfile);
if fidout1==-1
    error(helpword);
end
for c1=1:(M+1)
    fprintf(fidout1,'%6.4f %6.4f %6.4f\n',T(c1),sk(c1),standard(c1));
end
fclose(fidout1);
fidout2=fopen('illuminate.txt','w');
if fidout2==-1
    error('Can''t open the file----illuminate.txt');
end
fprintf(fidout2,'Sequence length:%6.0f\n',N);
fprintf(fidout2,'Value of M:%6.0f\n',M);
fprintf(fidout2,'Illumination of every rows:\n');
fprintf(fidout2,'    1--Periods    2--The smoothed power spectra    3--The standard noice spectra.');
fclose(fidout2);
disp('All done!');