KPowermeans算法进行毫米波通信中多径分量聚类，结合CH,DB索引来找最佳簇数

function [CH_Value,DB_Value] = CombinedValidate(Indexdata,ClusterCenter,Data) %COMBINEDVALIDATE Summary of this function goes here % Detailed explanation goes here % % % 输入参数 ---Indexdata 此参数是一次kpowermeans算法的运行结果，每个数据的类索引(num*2) % ---ClusterCenter kpowermeans算法输出的类中心 % ---Data kpowermeans算法利用的原始数据 % 输出参数 --- CH_Value CH算法算出来的聚类结果的有效性度量值 % ---DB_Value DB计算出来的类结果的有效性度量值 %获取数据的个数number和簇的数目k和参数的维度cols %求得每个簇内多径数目L_k(k*1) %获得原始数据data （AOA_a,AOA_e,AOD_a,AOD_e,DELAY,POWER）（num*6） %获得各个类的中心 [Number,cols]=size(Data); [k,~]=size(ClusterCenter); %如果不加~，那么输出k就是一个向量 L_k=zeros(k,1); data1=Data; center=ClusterCenter; centercopy=ClusterCenter; %计算距离时后必要的引入一些计算量 % MCD(delay)(ij)=5*|delay(i)-delay(j)|*delay(std)/delay(max)*delay(max) % MCD(AOA/AOD)(ij)=1/2* |()-()| min_delay=min(data1(:,5)); max_delay=max(data1(:,5)); std_delay=std(data1(:,5)); angle1_aoa=[0,0,0]; % 第i和j个数据点的角度度量 angle2_aoa=[0,0,0]; angle1_aod=[0,0,0]; angle2_aod=[0,0,0]; MCD_delay=0; % 第i和j个数据点的时延度量 MCD_AOA=0; MCD_AOD=0; MCD_POWER=0; %求整个数据集的中心 power_total=sum(data1(:,6)); c=zeros(1,6); c(1)=dot(data1(:,6),data1(:,1))/power_total; c(2)=dot(data1(:,6),data1(:,2))/power_total; c(3)=dot(data1(:,6),data1(:,3))/power_total; c(4)=dot(data1(:,6),data1(:,4))/power_total; c(5)=dot(data1(:,6),data1(:,5))/power_total; c(6)=power_total/Number; %计算每个类里面多径的数目 for i=1:Number L_k(Indexdata(i,2),1)=L_k(Indexdata(i,2),1)+1; end %利用CH索引来计算k分类下的索引value %利用CH索引来计算k分类下的索引value %利用CH索引来计算k分类下的索引value %类间计算 MCD_IN=0; MCD_OUT=0; for j=1:k MCD_delay=15*abs(centercopy(j,5)-c(5))*std_delay/(max_delay-min_delay)*(max_delay-min_delay); angle1_aoa=[sin(centercopy(j,1))*cos(centercopy(j,2)),sin(centercopy(j,1))*sin(centercopy(j,2)),cos(centercopy(j,1))]; angle2_aoa=[sin(c(1))*cos(c(2)),sin(c(1))*sin(c(2)),cos(c(1))]; angle1_aod=[sin(centercopy(j,3))*cos(centercopy(j,4)),sin(centercopy(j,3))*sin(centercopy(j,4)),cos(centercopy(j,3))]; angle2_aod=[sin(c(3))*cos(c(4)),sin(c(3))*sin(c(4)),cos(c(3))]; MCD_AOA=0.5*0.5*sum((angle1_aoa-angle2_aoa).^2); MCD_AOD=0.5*0.5*sum((angle1_aod-angle2_aod).^2); %MCD_POWER=abs(centercopy(j,6)-data1(i,6)); dist=L_k(j,1)*(MCD_AOA+MCD_AOD+MCD_delay*MCD_delay); MCD_OUT=MCD_OUT+dist; end %类内计算 for i=1:Number MCD_delay=15*abs(centercopy(Indexdata(i,2),5)-data1(i,5))*std_delay/(max_delay-min_delay)*(max_delay-min_delay); angle1_aoa=[sin(centercopy(Indexdata(i,2),1))*cos(centercopy(Indexdata(i,2),2)),sin(centercopy(Indexdata(i,2),1))*sin(centercopy(Indexdata(i,2),2)),cos(centercopy(Indexdata(i,2),1))]; angle2_aoa=[sin(data1(i,1))*cos(data1(i,2)),sin(data1(i,1))*sin(data1(i,2)),cos(data1(i,1))]; angle1_aod=[sin(centercopy(Indexdata(i,2),3))*cos(centercopy(Indexdata(i,2),4)),sin(centercopy(Indexdata(i,2),3))*sin(centercopy(Indexdata(i,2),4)),cos(centercopy(Indexdata(i,2),3))]; angle2_aod=[sin(data1(i,3))*cos(data1(i,4)),sin(data1(i,3))*sin(data1(i,4)),cos(data1(i,3))]; MCD_AOA=0.5*0.5*sum((angle1_aoa-angle2_aoa).^2); MCD_AOD=0.5*0.5*sum((angle1_aod-angle2_aod).^2); %MCD_POWER=abs(centercopy(j,6)-data1(i,6)); MCD_IN=MCD_IN+(MCD_AOA+MCD_AOD+MCD_delay*MCD_delay); end %计算CH value CH_Value=MCD_OUT*(Number-k)/(k-1)/MCD_IN; disp('CH_Value'); disp(CH_Value); %利用DB索引来计算k分类下的索引value %利用DB索引来计算k分类下的索引value %利用DB索引来计算k分类下的索引value Sk=zeros(k,1); for i=1:Number MCD_delay=15*abs(centercopy(Indexdata(i,2),5)-data1(i,5))*std_delay/(max_delay-min_delay)*(max_delay-min_delay); angle1_aoa=[sin(centercopy(Indexdata(i,2),1))*cos(centercopy(Indexdata(i,2),2)),sin(centercopy(Indexdata(i,2),1))*sin(centercopy(Indexdata(i,2),2)),cos(centercopy(Indexdata(i,2),1))]; angle2_aoa=[sin(data1(i,1))*cos(data1(i,2)),sin(data1(i,1))*sin(data1(i,2)),cos(data1(i,1))]; angle1_aod=[sin(centercopy(Indexdata(i,2),3))*cos(centercopy(Indexdata(i,2),4)),sin(centercopy(Indexdata(i,2),3))*sin(centercopy(Indexdata(i,2),4)),cos(centercopy(Indexdata(i,2),3))]; angle2_aod=[sin(data1(i,3))*cos(data1(i,4)),sin(data1(i,3))*sin(data1(i,4)),cos(data1(i,3))]; MCD_AOA=0.5*0.5*sum((angle1_aoa-angle2_aoa).^2); MCD_AOD=0.5*0.5*sum((angle1_aod-angle2_aod).^2); Sk(Indexdata(i,2),1)=Sk(Indexdata(i,2),1)+sqrt(MCD_AOA+MCD_AOD+MCD_delay*MCD_delay); end for i=1:k Sk(k,1)=Sk(k,1)/L_k(k,1); end d_ij=zeros(k,k); for i=1:k for j=1:k MCD_delay=15*abs(centercopy(i,5)-centercopy(j,5))*std_delay/(max_delay-min_delay)*(max_delay-min_delay); angle1_aoa=[sin(centercopy(i,1))*cos(centercopy(i,2)),sin(centercopy(i,1))*sin(centercopy(i,2)),cos(centercopy(i,1))]; angle2_aoa=[sin(centercopy(j,1))*cos(centercopy(j,2)),sin(centercopy(j,1))*sin(centercopy(j,2)),cos(centercopy(j,1))]; angle1_aod=[sin(centercopy(i,3))*cos(centercopy(i,4)),sin(centercopy(i,3))*sin(centercopy(i,4)),cos(centercopy(i,3))]; angle2_aod=[sin(centercopy(j,3))*cos(centercopy(j,4)),sin(centercopy(j,3))*sin(centercopy(j,4)),cos(centercopy(j,3))]; MCD_AOA=0.5*0.5*sum((angle1_aoa-angle2_aoa).^2); MCD_AOD=0.5*0.5*sum((angle1_aod-angle2_aod).^2); d_ij(i,j)=sqrt(MCD_AOA+MCD_AOD+MCD_delay*MCD_delay); end end R=zeros(k,1); for i=1:k R(i,1)=-inf; d_ij_copy= d_ij; %重新复制数组，否则我们在不断改变实际的数组值； for j=1:k if j==i data=-5000; if data> R(i,1) R(i,1)=data; end else data=(Sk(i,1)+Sk(j,1))/d_ij_copy(i,j); if data> R(i,1) R(i,1)=data; end end end end % for i=1:k % datamatix=Sk(i,1)*ones(1,k-1); % data1matrix=Sk; % data1matrix(i)=[]; % data2=d_ij; % data2(:,i)=[]; % R(i,1)=max((datamatix+data1matrix(:,1)')/data2(i,:)); % end %%最后计算DB值 DB_Value=sum(R(:,1))/k; disp('DB_Value'); disp(DB_Value); end