cdmid.m

来自「signal and system answer」· M 代码 · 共 77 行

clear
load cdmid
nihex2=zeros(111,8);
nihey2=ones(111,8)*1000;
for i=1:111
    k=1;
    for j=1:45
        if cdfen2nihe(i,j)~=1000&xiabiao2(i,j)~=0
            nihex2(i,k)=xiabiao2(i,j);
            nihey2(i,k)=cdfen2nihe(i,j);
            k=k+1;
        end
    end
end
cdfen2nihe=0.2*cdfen2nihe;
nihey2=0.2*nihey2;
l=ones(1,111);
for i=1:111
    k=1;
    while nihex2(i,k)~=0
       k=k+1;
    end
    k=k-1;
    x2=nihex2(i,1:k);
    while nihey2(i,l(i))~=1000*0.2
       l(i)=l(i)+1;
    end
    l(i)=l(i)-1;
    y2=nihey2(i,1:l(i));
   
    pcd2(i,:)=polyfit(x2,y2,2);
end
for i=1:111
    for j=1:45
        if cdfen2nihe(i,j)==1000*0.2;
            cdfen2nihe(i,j)=polyval(pcd2(i,:),j);
        else
            l(i)=l(i)-1;
        end
          if cdfen2nihe(i,j)<=0|l(i)==0
            cdfen2nihe(i,j)=0;
        end
    end
end
bilv2cd=zeros(111,45);
for j=1:111
    k=1;
    while k<=45&cdfen2nihe(j,k)~=0
        k=k+1;
    end
    k=k-1;
    for i=1:k
        if i==1
            bilv2cd(j,1)=(-3*cdfen2nihe(j,1)+4*cdfen2nihe(j,2)-cdfen2nihe(j,3))/2;
        elseif i~=k
            bilv2cd(j,i)=(cdfen2nihe(j,i+1)-cdfen2nihe(j,i-1))/2;
        elseif k>2
            bilv2cd(j,k)=(cdfen2nihe(j,k-2)-4*cdfen2nihe(j,k-1)+3*cdfen2nihe(j,k))/2;
        end
    end
end
numx2=zeros(1,45);numy2=zeros(1,45);
for i=1:111
    for j=1:45
        if bilv2cd(i,j)~=0
            numx2(j)=numx2(j)+1;
        end
        if cdfen2nihe(i,j)~=0
            numy2(j)=numy2(j)+1;
        end
    end
end
sc=sum(cdfen2nihe);
cdfen2nihemean=sc(1:44)./numy2(1:44);
month=1:44;
figure
plot(month,cdfen2nihemean,'r.')