PT_extaction.zip_matlab峰值提取_光谱峰值提取_谱峰提取

clc; clear all; %846 %guangpu=[3120 3116 3124 3134 3156 3146 3148 3177 3140 3118 3118 3119 3130 3134 3130 3157 3148 3150 3162 3180 3168 3126 3117 3126 3130 3144 3150 3156 3146 3158 3166 3166 3178 3136 3130 3118 3127 3138 3132 3160 3154 3164 3176 3168 3162 3140 3126 3134 3142 3134 3152 3154 3168 3170 3164 3180 3174 3136 3124 3118 3130 3122 3144 3149 3154 3158 3168 3180 3168 3130 3126 3118 3126 3126 3138 3140 3148 3166 3156 3164 3182 3160 3138 3124 3116 3124 3131 3142 3128 3138 3177 3156 3154 3146 3136 3119 3118 3123 3130 3137 3145 3166 3176 3168 3183 3172 3166 3138 3120 3130 3130 3124 3138 3138 3142 3150 3144 3176 3158 3136 3117 3118 3123 3127 3132 3124 3134 3159 3152 3164 3177 3154 3144 3138 3128 3122 3132 3136 3150 3158 3175 3166 3174 3185 3180 3140 3132 3118]; d=load('D:\DS\MATLAB\gooddata\109.txt'); %guangpu=d(839,151:497); %guangpu=d(844,250:397); %844 %guangpu=[3138 3136 3150 3158 3172 3164 3182 3170 3148 3140 3132 3140 3146 3150 3174 3158 3166 3168 3186 3188 3142 3128 3136 3134 3136 3146 3154 3172 3164 3168 3178 3186 3156 3130 3152 3148 3158 3156 3168 3177 3177 3170 3190 3178 3162 3150 3131 3154 3164 3162 3178 3186 3184 3194 3194 3174 3156 3146 3133 3144 3148 3150 3158 3172 3186 3190 3198 3182 3150 3150 3136 3146 3156 3146 3158 3168 3174 3162 3182 3180 3166 3156 3132 3128 3140 3152 3142 3132 3156 3154 3158 3158 3170 3144 3118 3123 3124 3119 3123 3128 3150 3164 3189 3190 3189 3182 3146 3142 3142 3138 3130 3144 3150 3144 3164 3164 3184 3176 3150 3128 3126 3132 3126 3140 3142 3156 3168 3152 3186 3184 3188 3156 3132 3134 3142 3150 3158 3158 3176 3176 3180 3178 3190 3192 3168 3148 3120]; wavenum=1024;%波长数 V_src=[];%变异系数 V_PT=[]; M_src=[];%均值 M_PT=[]; S_src=[];%标准差 S_PT=[]; DS_DATA=[];%存放动态光谱 DS_DATA_PT=[]; %设置进度条，查看提取进度 handleofwaitbar=waitbar(0); for w_n=1:wavenum guangpu_src=d(w_n,151:397); guangpuPT=[]; gzzb=[]; %谷点坐标 fzzb=[]; %峰值坐标 zqcd=[]; %周期长度 gz=0;%谷值数量 fz=0;%峰值数量 %求取理想波形的波峰和波谷值的坐标，用来提取前沿后沿模板 %最大程度的滤波，便于极值点的挑选 %选择参数 if guangpu_src(1)<500 %可见光 Wp=0.2; %通带截止频率，为归一化的频率，在0~1之间，1对应抽样频率的一半,此处即50Hz Ws=0.4; %阻带截止频率，此处为10Hz Rp=0.01; %通带的衰减dB Rs=50; %阻带的衰减dB else %近红外 Wp=0.3; Ws=0.5; Rp=0.01; Rs=50; end %滤波 [N,Wn]=buttord(Wp,Ws,Rp,Rs); %其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。 [B,A]=butter(N,Wn); %Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。 guangpu=filtfilt(B,A,guangpu_src); %零相位数字滤波器，可使起始和结束阶段的暂态过度过程最短 %寻找模板 %窗口滑动寻找谷值点 Ns=length(guangpu); %总的数据点个数 valley=0; temp=0; for k=8:(Ns-7) if k==8 %when then start point is not a valley,this if-loop is unnecessary [C,I]=min(guangpu(k-7:k+7)); if I<=k gzzb=[gzzb I]; gz=gz+1; if gz>1 temp=gzzb(gz)-gzzb(gz-1); zqcd=[zqcd temp]; end end elseif k==(Ns-7) %when then start point is not a valley,this if-loop is unnecessary [C,I]=min(guangpu(k-7:k+7)); I=(k-7)+(I-1); % I为索引值 将其转变为真实值 if I>=k&&(I-gzzb(end)>5) gzzb=[gzzb I]; gz=gz+1; if gz>1 temp=gzzb(gz)-gzzb(gz-1); zqcd=[zqcd temp]; end end else if guangpu(k)==min(guangpu(k-7:k+7)) if gz==0||(k-gzzb(end)>5) % valley=find(guangpu(k-8:k+8)==min(guangpu(k-8:k+8))) gzzb=[gzzb k]; gz=gz+1; if gz>1 temp=gzzb(gz)-gzzb(gz-1); zqcd=[zqcd temp]; end end end end end %窗口滑动寻找峰值点 for k=(gzzb(1)+7):(gzzb(end)-7) %如果在起始处 if k==(gzzb(1)+7) [C,I]=max(guangpu(k-7:k+7)); I=(k-7)+(I-1); % I为索引值 将其转变为真实值 if I<k fzzb=[fzzb,I]; fz=fz+1; end end %如果滑动到(gzzb(end)-7) if k==(gzzb(end)-7); [C,I]=max(guangpu(k-7:k+7)); I=(k-7)+(I-1); % I为索引值 将其转变为真实值 if I>k&&(k-fzzb(end)>5) fzzb=[fzzb,I]; fz=fz+1; end end %判断是否为峰值 if guangpu(k)==max(guangpu(k-7:k+7)) if fz==0||(k-fzzb(end)>5) % valley=find(guangpu(k-8:k+8)==min(guangpu(k-8:k+8))) fzzb=[fzzb k]; fz=fz+1; end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %分别对各个周期(谷值之间)进行峰值变换 for j=1:(length(gzzb)-1) %每次清空 zeng=[]; jian=[]; guangpuPT(gzzb(j))=guangpu(gzzb(j)); if j==(length(gzzb)-1) guangpuPT(gzzb(end))=guangpu(gzzb(end)); end for n=gzzb(j):(gzzb(j+1)-1) delta=guangpu(n+1)-guangpu(n); if delta>=0 zeng=[zeng delta]; else jian=[jian delta]; end end zeng=[zeng jian]; leijia=0; count=1; for n=gzzb(j):(gzzb(j+1)-1) leijia=leijia+zeng(count); count=count+1; guangpuPT(n+1)=guangpu(gzzb(j))+leijia; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %PT变换后的数据的谷值点不变 不需要再次搜索 gzzb[] %求原始光谱中幅值存入vpp vpp_src=[]; if fz<gz for i=1:fz temp=guangpu(fzzb(i))-guangpu(gzzb(i)); vpp_src=[vpp_src temp]; temp=guangpu(fzzb(i))-guangpu(gzzb(i+1)); vpp_src=[vpp_src temp]; end end if fz>gz for i=2:gz temp=guangpu(fzzb(i))-guangpu(gzzb(i-1)); vpp_src=[vpp_src temp]; temp=guangpu(fzzb(i))-guangpu(gzzb(i)); vpp_src=[vpp_src temp]; end end if fz==gz if fzzb(1)<gzzb(1) for i=2:gz temp=guangpu(fzzb(i))-guangpu(gzzb(i)); vpp_src=[vpp_src temp]; temp=guangpu(fzzb(i))-guangpu(gzzb(i-1)); end elseif fzzb(1)>gzzb(1) for i=1:gz-1 temp=guangpu(fzzb(i))-guangpu(gzzb(i)); vpp_src=[vpp_src temp]; temp=guangpu(fzzb(i))-guangpu(gzzb(i+1)); end end end %窗口滑动寻找峰值点 fzzbPT=[]; %PT变换后的峰值坐标 fzPT=0; %PT变换后的峰值数量 %窗口滑动寻找峰值点 for k=(gzzb(1)+6):(gzzb(end)-6) %如果在起始处 if k==(gzzb(1)+6) [C,I]=max(guangpuPT(k-6:k+6)); I=(k-6)+(I-1); % I为索引值 将其转变为真实值 if I<k fzzbPT=[fzzbPT,I]; fzPT=fzPT+1; end end %如果滑动到(gzzb(end)-7) if k==(gzzb(end)-6); [C,I]=max(guangpuPT(k-6:k+6)); I=(k-6)+(I-1); % I为索引值 将其转变为真实值 if I>k&&(k-fzzbPT(end)>5) fzzbPT=[fzzbPT,I]; fzPT=fzPT+1; end end %判断是否为峰值 if guangpuPT(k)==max(guangpuPT(k-6:k+6)) if fzPT==0||(k-fzzbPT(end)>5) % valley=find(guangpu(k-8:k+8)==min(guangpu(k-8:k+8))) fzzbPT=[fzzbPT k]; fzPT=fzPT+1; end end end %求PT变换后光谱中幅值存入vpp_PT vpp_PT=[]; if fzPT<gz for i=1:fzPT temp=guangpuPT(fzzbPT(i))-guangpuPT(gzzb(i)); vpp_PT=[vpp_PT temp]; temp=guangpuPT(fzzbPT(i))-guangpuPT(gzzb(i+1)); vpp_PT=[vp