readdata.m

本次的算法是相位周期展开聚类

% deal data from specified data file
clear; 
mird=0.2;
data_fname1='111.bin'; % 这里是文件名
file_id1=fopen(data_fname1,'rb'); 
data_fname2='222.bin'; % 这里是文件名
file_id2=fopen(data_fname2,'rb'); 
s1=[zeros(1,128)];
s3=[zeros(1,128)];
%bmaxid=fopen('1111.bin','a');
%dmaxid=fopen('2222.bin','a');
while feof(file_id1)~=1&feof(file_id2)~=1 
    % 数据类型是int16，每次读入256个数
    % raw_array每次都是256x1的矩阵，ele_count为读入的数的个数（正常情况下应为256）
    [row_array1,ele_count]=fread(file_id1,128,'double'); 
    [row_array2,el_count]=fread(file_id2,128,'double'); 
    if ele_count <128 % elecount < 512代表数据不够，已经到了文件的结尾
        break; 
    else 
        % 将256x1的row_array转置一下，变为1x256的矩阵
        % s1=[zeros(128)];
        %s2=[s1,s(1:128)];
        % y=gccphat(t1,t2);
        % s1=s(129:256);
        row_array1=row_array1'; 
        s2=[s1,row_array1(1:128)];
        s1=row_array1(1:128);
    %[row_array2, el_count]=fread(file_id2,128, 'double'); 
    if el_count <128 % elcount < 512代表数据不够，已经到了文件的结尾
        break; 
    else 
        row_array2=row_array2'; 
        s4=[s3,row_array2(1:128)];
        s3=row_array2(1:128);
        nfft=256;
        f1=fft(s2,nfft);
        %fx1=f1(1:256);
        f2=fft(s4,nfft);
        %fx2=f2(1:256);
        fz2=conj(f2);
        %nfft=256;
        %[cxy,f]=cpsd[s2,s4,nfft,hanning(nfft)];
        Pxy=f1.*fz2;
        %Pxy=Pxy./sqrt(Pxy.*conj(Pxy));
        %Pxy=ifft(R);
        delta=angle(Pxy);
        da=1:256;
        f=1/16000;
        w=da/(256*f);        
        %dl=Ap./w;
        %M=pin(Ap,dl);
        sun=zeros(128,6);
        %sun=[];
        M=[];
        for i=1:128
           dsum=dlay(delta(i),mird,w(i));
           dsum=dsum*16000;
           sun(i,:)=dsum;
           sun=sun;
           %sun=[sun,dsum];
           %sun=abs(sun);
        end
        for m=1:128
            for j=1:6
            if (sun(m,j)~=0)
                M=[M,sun(m,j)];
            end
            end
        end
         %sun1=sort(M,'ascend');
         %b=unique(sun1);
         %c=zeros(size(b));
         %for n=1:length(b);
             %c(n)=length(find(sun1==b(n)));
         %end
         b=zeros(1,ceil(max(M)+ceil(abs(min(M)))));
c=[];
for n=ceil(min(M)):ceil(max(M))
c=[c,n];
end   
for m=1:length(M)
    for n=1:length(c)
        if (M(m)>=c(n)&M(m)<c(n)+1)
            b(n)=b(n)+1;
        end
    end
end
[dmax,imax]=max(b);
s=zeros(dmax,6);
%for i=1:dmax
ie=1;
wp=[];
   for m=1:128
      for n=1:6
        if(sun(m,n)>=c(imax)&sun(m,n)<(c(imax)+1)&sun(m,n)~=0)
         s(ie,:)=sun(m,:);
         wp=[wp,62.5*m];
         ie=ie+1;
        %end
        end          
    end
end
        end
    end
end       
%fclose(bmaxid);
%fclose(dmaxid);
fclose(file_id1);
fclose(file_id2);
% 关闭文件
clear data_fname file_id fid ele_count i m n row_array ans;