从fit文件中提取lick指数的matlab程序

lick指数的定义及公式，可以参考以下网址：

http://d.wanfangdata.com.cn/periodical/dlxyxb201406010

本程序实现从fit文件中读取lick指数

clear;
clc;
tic;readDir = 'E://cas/b6';    %direction of fitsreadPath = [readDir '\*.fit'];readList = dir(readPath);[m1, n1] = size(readList);
dim_num=3522;                                                       %total dimensions
x=0; k1=1;%0-5k2=1;%5-10k3=1;%10-15k4=1;%15-20k5=1;%20-25k6=1;%25-30k7=1;%30-40k8=1;%40-60k9=1;%60-80k10=1;k11=1;k12=1;for i1 = 1:m1                                                        %templetefilename =[readDir '\'  readList(i1, 1).name];    %name of file%splate1=fitsread(picName);                                      %read data%spec=splate1(1,:);                                              %first line% filename='E:\cas\b6\spSpec-51691-0342-608.fit'info = fitsinfo(filename);
l = info.PrimaryData.Keywords;j = 1;
s = '';
while(j <= length(l))
s = strvcat(s,char(l(j,1)));
j = j + 1;
end%步长
SN_G   = strmatch('SN_G',s,'exact');
sn(i1)=l(SN_G,2);   %提取信噪比,sn为原始的信噪比数组SN=cell2mat(sn(i1));   if(isnumeric(SN))Sn=SN;                                                         %提取数值型的信噪比
elseSn=-1;
end
sn2(1,i1)=Sn(1,1);
%sn2为最终可用的信噪比,为一行m1列的向量，m1为fit文件个数%这样，针对一个特定fit文件，我们得到了他的信噪比NAXIS1 = strmatch('NAXIS1',s,'exact');    %fits文件的维度，注意，很多fits文件的维度是不同的
NAXIS1=l(NAXIS1,2);
NAXIS1=cell2mat(NAXIS1);
COEFF0   = strmatch('COEFF0',s,'exact');  %起始波长,注意1，不同文件起始波长也往往是不一样的
COEFF0  =l(COEFF0  ,2);                  %注意2，设该值为x,则10进制的起始波长为10^x
COEFF0  =cell2mat(COEFF0  );
COEFF1   = strmatch('COEFF1',s,'exact');  %步长，注意，也是取了log,设该值为x,则10进制的步长为10^x
COEFF1  =l(COEFF1  ,2);                  %第i个点的波长为，10^(COEFF0 + COEFF1*i),  具体参考见：http://classic.sdss.org/dr2/products/spectra/read_spSpec.html
COEFF1  =cell2mat(COEFF1);
%--解析头文件里内容--endsplate1=fitsread(filename);              %一般对于dr8及以前的fits文件，用默认的fitsread()函数即可
spec=splate1(1,:);                       %第一行是真正的流量，如果读第x行，就是 splate1(x,:);  %---波长数组--begin
wave=ones(1,NAXIS1);                     %生成一个长度为NAXIS1的一维数组，用于存储波长
for i=1:NAXIS1   wave(i)=i-1;  end;      %matlab数组的下标从1开始，而不是0；数组的值为0,1,2,3,4........
logwavelength = COEFF0 + wave * COEFF1;   %开始处理波长数组
for i=1:NAXIS1   wave(i)=10^logwavelength(i);  end;  % 转变为10进制的波长，第i个点的10进制波长就是wave(i)
%---波长数组--endindex_bind_L = [4222.250,4281.375,4369.125,4452.125,4514.250,4634.000,4847.875,4977.750,5160.125,5245.650, 5312.125, 5387.500, 5696.625, 5776.625, 5876.875, 4084.750,4321.000,4092.250, 4332.500,4142.125,4142.125,5069.125,5154.125,5936.625,6189.625];
index_bind_R = [4234.750,4316.375,4420.375,4474.625,4559.250,4720.250,4876.625,5054.000,5192.625,5285.650,5352.125,5415.000,5720.375,5796.625,5909.375,4123.500,4364.750,4113.500,4353.500,4177.125,4177.125,5134.125,5196.625,5994.125,6272.125];
blue_L = [4211.000,4266.375,4359.125,4445.875,4504.250,4611.500,4827.875,4946.500,5142.625,5233.150,5304.625,5376.250,5672.875,5765.375,5860.625,4042.850,4284.750,4058.500,4284.750,4080.125,4083.875,4895.125,4895.125,5816.625,6066.625];
blue_R = [4219.750,4282.625,4370.375,4454.625,4514.250,4630.250,4847.875,4977.750,5161.375,5248.150,5315.875,5387.500,5696.625,5775.375,5875.625,4081.000,4321.000,4089.750,4321.000,4117.625,4096.375,4957.625,4957.625,5849.125,6141.625];
red_L = [4241.000,4318.875,4442.875,4477.125,4560.500,4742.750,4876.625,5054.000,5191.375,5285.650,5353.375,5415.000,5722.875,5797.875,5922.125,4129.750,4368.500,4116.000,4356.000,4244.125,4244.125,5301.125,5301.125,6038.625,6372.625];
red_R = [4251.000,4335.125,4455.375,4492.125,4579.250,4756.500,4891.625,5065.250,5206.375,5318.150,5363.375,5425.000,5736.625,5811.625,5948.125,4162.250,4421.000,4138.500,4386.000,4284.125,4284.125,5366.125,5366.125,6103.625,6415.125];%25个区间的端点
index_L0 = get_index_function(wave,index_bind_L);
index_R0 = get_index_function(wave,index_bind_R);
blue_L0 = get_index_function(wave,blue_L);
blue_R0 = get_index_function(wave,blue_R);
red_L0 = get_index_function(wave,red_L);
red_R0 = get_index_function(wave,red_R);licks = zeros(1,25);
for i=1:19                                   %在一个带中
wave_Mid=wave(index_L0(i):index_R0(i));
wave_Blue=wave(blue_L0(i):blue_R0(i));
wave_Red=wave(red_L0(i):red_R0(i)); Mid_flux=spec(index_L0(i):index_R0(i));
Blue_flux=spec(blue_L0(i):blue_R0(i));
Red_flux=spec(red_L0(i):red_R0(i));lamda1_b=min(wave_Blue);
lamda2_b=max(wave_Blue);
lamda1_r=min(wave_Red);
lamda2_r=max(wave_Red);
lamda1_mid=min(wave_Mid);
lamda2_mid=max(wave_Mid);
mean_flux_b (1,i)= (lamda1_b+lamda2_b)/2;
mean_flux_b(2,i) = trapz(wave_Blue,Blue_flux)/(lamda2_b-lamda1_b);
mean_flux_r(1,i)=(lamda1_r+lamda2_r)/2;
mean_flux_r(2,i)=trapz(wave_Red,Red_flux)/(lamda2_r-lamda1_r);%flux_mid是一个向量，存储中间带积分值
% red_flux(i) = get_flux (wave_Red(i),wave,spec);
% blue_flux(i) = get_flux(wave_Blue,wave,spec);  %由图中点对应的离散的横坐标（波长）temp1 =mean_flux_r(1,i)-mean_flux_b(1,i);%x2-x1temp2 = mean_flux_r(2,i)-mean_flux_b(2,i);%y2-y1tempk = temp2/temp1;FClamda = mean_flux_b(2,i)+tempk*(wave_Mid-mean_flux_b(1,i));if i<=19  licks(1,i) =real(double(trapz(wave_Mid,(1-Mid_flux./FClamda))));elselicks(1,i) =real(double(-2.5*log(trapz(wave_Mid,(Mid_flux./FClamda)/(lamda2_mid-lamda1_mid)))));endend%这样，针对一个fit文件，我们得到了它的25维lick数组licksendtoc;

从fit文件中提取lick指数的matlab程序相关推荐

ML之MLiR：利用多元线性回归法，从大量数据(csv文件)中提取五个因变量(输入运输任务总里程数、运输次数、三种不同的车型，预测需要花费的小时数)来预测一个自变量
ML之MLiR:利用多元线性回归法,从大量数据(csv文件)中提取五个因变量(输入运输任务总里程数.运输次数.三种不同的车型,预测需要花费的小时数)来预测一个自变量输出结果代码设计 from nu ...
[SimplePlayer] 4. 从视频文件中提取音频
提取音频,具体点来说就是提取音频帧.提取方法与从视频文件中提取图像的方法基本一样,这里仅列出其中的不同点: 1. 由于目的提取音频,因此在demux的时候需要指定的是提取audio stream Au ...
python提取文件指定列_如何从csv文件中提取特定列并使用python绘图
我有一个csv文件,其中包含以下几行数据:# Vertex X Y Z K_I K_II K_III J 0 2.100000e+00 2.000000e+00 -1.000000e-04 0.000 ...
gnuradio上怎么使用python文件_使用Python从PDF文件中提取数据
前言数据是数据科学中任何分析的关键,大多数分析中最常用的数据集类型是存储在逗号分隔值(csv)表中的干净数据.然而,由于可移植文档格式(pdf)文件是最常用的文件格式之一,因此每个数据科学家都应该了 ...
ffmpeg-从mp4、flv、ts文件中提取264视频流数据
ffmpeg-从mp4.flv.ts文件中提取264视频流数据 main.c #include <stdio.h> #include <libavutil/log.h> #in ...
ffmpeg-从flv文件中提取AAC音频数据保存为文件
AAC ADTS格式协议: 从flv文件中提取AAC音频数据保存为文件. 如果需要详细了解AAC ADTS格式,可以查询文档. 原文件: 提取aac文件: main.c #include <st ...
java 取pdf 文本域_java – 使用iText从pdf文件中提取文本列
我需要使用iText从pdf文件中提取文本. 问题是:一些pdf文件包含2列,当我提取文本时,我得到一个文本文件,其中列被合并为结果(即同一行中两列的文本) 这是代码: public class pd ...
一行命令从 APK 文件中提取 Endpoint 及 URL
做IoT的人免不了要接触Android,接触Android的人又免不了要研究别人的App应用. Diggy,一款能够从 apk 文件中提取 endpoint 及 URL 的工具,只要一行命令就可以帮大 ...
使用Python从PDF文件中提取数据
前言数据是数据科学中任何分析的关键,大多数分析中最常用的数据集类型是存储在逗号分隔值(csv)表中的干净数据.然而,由于可移植文档格式(pdf)文件是最常用的文件格式之一,因此每个数据科学家都应该了 ...
matlab出如何从fig中获取数据,如何从MATLAB的fig文件中提取原始数据？
如何从MATLAB的fig文件中提取原始数据? mip版关注:171 答案:3 悬赏:70 解决时间 2021-02-23 07:29 已解决 2021-02-23 02:41 如何从MATL ...

从fit文件中提取lick指数的matlab程序

从fit文件中提取lick指数的matlab程序相关推荐

最新文章

热门文章