clear;
P1=load('200540-1.txt')'
PP=P1(351:450);
PS=interp(PP,10);%10內插
P=medfilt1(PS,9);%中值濾波
np=min(P); %信號最小值
mp=max(P); %信號最大值
dp=np+(mp-np)./4; %信號下四分之一
up=np+2.*(mp-np)./3;%信號上三分之二
%P=medfilt1(Pp,3);
N=[];
for n=2:999;
    pc1=P(n)-P(n-1);
    pc2=P(n)-P(n+1);
    pcc=pc1.*pc2;
    if pcc>=0;
        N=[N,n];%把符合條件的峰谷像素值放進數組N中
    end  
end
s=length(N);%計算數組N的長度
%N1=N(2:s);% 尋訪從第2個元素到最后一個元素，并將其放入數組N1中
px=0; %存放不獨立的峰谷位置標志
for i=2:s
    xp=N(i)-N(i-1);
    pp=P(N(i))-P(N(i-1));
    if (xp==1)&(abs(pp)==0) %該峰谷位置比前一個大一個像素，且灰度值相等，則該峰（谷）位置不獨立
        px(i)=1; %不獨立的峰谷位置標志為1
    elseif (abs(pp)<=2)&((P(N(i))<dp)|(P(N(i))>up))%灰度值大于up或小于dp的位置，應為峰谷區，此時灰度值有波動或峰谷位置不相鄰也可能是不獨立的峰谷;(灰度值處于up和dp之間的波動應視為峰谷波動）
         px(i)=1;
    else
        px(i)=0;
    end
end
for j=1:s-1
    dpx(j)=px(j)-px(j+1);
end
apx=find(dpx==-1);
bpx=find(dpx==1);
lapx=length(apx);
lbpx=length(bpx);
lpx=min(lapx,lbpx);
N2=[];
tp=0;
for ii=1:s-1
    if (px(ii)==0)&(dpx(ii)==0)
       N2=[N2,N(ii)];
    elseif (px(ii)==0)&(dpx(ii)==-1)&(tp<=lpx)
       uua=N(ii);%不獨立峰谷位置的起始
    elseif (px(ii)==1)&(dpx(ii)==1)&(tp<=lpx)
       uub=N(ii);
       tp=tp+1; %不獨立峰谷位置的結束
       N0=fix((uua+uub)./2);
       N2=[N2,N0];
    end
end

lg=length(N2);%N2存儲的峰谷位置像素所組成的數組的長度

%下面找各峰谷中心位置坐標
w=P(N2);
M=[];
MM=[];
for i=2:lg-1
    T1=(w(i)+w(i-1))./2;%計算左閾值
    T2=(w(i)+w(i+1))./2;%計算右閾值
    t1=N2(i-1)+1;
    tm1=abs(P(t1)-T1);
    a=t1; %左閾值坐標
    while t1<N2(i)
          if abs(P(t1)-T1)<tm1;
              tm1=abs(P(t1)-T1);
              a=t1;
          else 
              t1=t1+1;
          end
    end 
    t2= N2(i)+1;
    tm2=abs(P(t2)-T2);
    b=t2; %右閾值坐標
    while t2<N2(i+1)
          if abs(P(t2)-T2)<tm2;
              tm2=abs(P(t2)-T2);
              b=t2;
          else 
              t2=t2+1;
          end
    end 
     S=0;
     T=0;
    for k=a:b
        S=S+k.*P(k);
        T=T+P(k);
    end
    m=S./T;%做加權平均的計算，得出中心位置
    MM=[MM,m];
    M=[M,fix(m)];%把中心位置存放入數組M中
end
%M=N2;
%plot(P),hold on
%plot(M,P(M),'*'),hold on 
%plot(N2,P(N2),'ro'),hold on
l=length(M);
ds=0;
% 下面計算條碼寬度
for i=2:1000
    ds(i)=P(i)-P(i-1); %信號求導
    %if ds(i)==0
        %ds(i)=1;
    %end
end
cds=find(ds==max(ds)); %導數最大值坐標
for ia=1:l-1
    if MM(ia)<cds(1)&cds(1)<MM(ia+1)
       as=ia;   %邊緣左邊峰（谷）位置
       bs=ia+1; %邊緣右邊谷（峰）位置
    end  
end  
S=0;
T=0;
for ks=M(as):M(bs)
        S=S+ks.*ds(ks);
        T=T+ds(ks);
end
e0=S./T; %找條碼的獨立邊緣e0
%計算各邊緣位置
E=[];%存儲各邊緣位置
E(as)=e0;
iq=as-1;
while iq>=1
     ae(iq)=2.*MM(iq+1);
     be(iq)=E(iq+1);
     E(iq)=ae(iq)-be(iq);
     iq=iq-1;
end
for ip=(as+1):l-1
    E(ip)=2.*MM(ip)-E(ip-1);
end 
FE=fix(E);
%計算各條碼寬度值
ll=l-1;%（l-1）個邊緣
wid=[];
if ds(FE(1))<0
     lc=0
     if rem(ll,2)==0
       la=ll./2;
     else 
       la=(ll-1)./2;
     end
else
     lc=1;
     if rem((ll-1),2)==0
       la=(ll-1)./2;
     else 
       la=(ll-2)./2;
     end
end
   
for lb=1:la
    lbb=2.*lb;
    wid(lb)=E(lbb+lc)-E(lbb-1+lc); %存儲各條碼的寬度值
end
MM=fix(M);
%subplot(2,1,1)
subplot(2,1,1),plot(P),hold on
plot(FE,P(FE),'ro')
subplot(2,1,2), plot(P),hold on
plot(M,P(M),'black*'),hold on
%plot(N2,P(N2),'r*'),hold on