drawingbandgap_3.m

一个用来进行频谱分析的程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%    本程序是参考了 J.Pendry 的程序选择了频谱密度函数作为确定本证频率的分析函数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
%%%%%%%%%%%%%%%%%%%%%%参数调整%%%%%%%%%%%%%%%%
LL1=20;   %这是峰的绝对高度判断参数
LL2=30;  %这是相邻点的峰相对差值判断参数 %%%这个参数的影响最大！！！！
LL3=2;   %这是合并靠得太紧的峰值的参数
fid=fopen('record_13_3.dat','r');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=fscanf(fid,'%f',1);
Dt=fscanf(fid, ' %f ',1);

c=fscanf(fid,'%f',1);

NTimeStep=fscanf(fid,'%f',1);
Nrecord=fscanf(fid,'%f',1);

Nofkx=fscanf(fid,'%f',1);
eigval=[];
 kk=[];
Ez=zeros(Nofkx,NTimeStep,Nrecord^2);
%%%% rea out the field Ez
for ii=1:Nofkx
   
    
    
    k(ii)=fscanf(fid,'%f',1);
    
    for jj=1:NTimeStep
           
         Ez(ii,jj,1:Nrecord^2)=fscanf(fid,'%f',Nrecord^2);
    
        end
end




for ii=1:Nofkx
    
    peaki1=[];
    sumspec=zeros(1,NTimeStep/2);
    Pxx=zeros(Nrecord^2,NTimeStep);
    spectrum=zeros(Nrecord^2,NTimeStep/2);
    for jj=1:Nrecord^2
        
        Pxx(jj,1:NTimeStep)=abs(fft(Ez(ii,:,jj)));
        
            spectrum(jj,1)=Pxx(jj,1)^2;
        
           for j1=2:NTimeStep/2
               
               spectrum(jj,j1)=Pxx(jj,j1)^2+Pxx(jj,NTimeStep+2-j1)^2;
               
           end
         
           for j1=1:NTimeStep/2
        
               spectrum(jj,j1)=spectrum(jj,j1)/NTimeStep^2;
               
           end
       end
       
       for j2=1:NTimeStep/2
           
         for jj=1:Nrecord^2
           
             sumspec(j2)=sumspec(j2)+spectrum(jj,j2);
            
         end
       end
  %   figure
  %  plot(sumspec);
  %  axis([0 60 0 max(sumspec)]);
     if sumspec(1)-sumspec(2)>max(sumspec)/LL2
                   
                     kk=[kk ii];
                     eigval=[eigval 0];
                 end
            
     
        for i1=2:NTimeStep/2-1
         
            if (sumspec(i1)-sumspec(i1-1)>max(sumspec)/LL2&sumspec(i1)-sumspec(i1+1)>max(sumspec)/LL2&sumspec(i1)>1/LL1*max(sumspec))
              %   if (sumspec(i1)-sumspec(i1-1)>max(sumspec)/LL2&sumspec(i1)-sumspec(i1+1)>max(sumspec)/LL2)
                
                 if length(peaki1)~=0
                             if i1-i0>LL3
               
                                     peaki1=[peaki1 i1];
                                      kk=[kk ii];
                                      w=(i1-1)/(Dt*NTimeStep)*a/(2*pi*c);
                                      eigval=[eigval w];
                             else  if sumspec(i1)>sumspec(i0)
                                         peaki1=[peaki1(1:length(peaki1)-1) i1];%????????????????
                                         % kk=[kk ii];?????????????????????
                                          w=(i1-1)/(Dt*NTimeStep)*a/(2*pi*c);%?????????????????????
                                          eigval=[eigval(1:length(eigval)-1) w]; % ??????????????????????             
                                     end
                             end
                 else
                         peaki1=[peaki1 i1]; 
                         kk=[kk ii];
                         w=(i1-1)/(Dt*NTimeStep)*a/(2*pi*c);
                         eigval=[eigval w];  
                  end
               i0=i1;
            end
        end
    %    cc=ii;
     

    end
        
      
figure(201)
plot(kk,eigval,'og')       
axis([1 16 0 0.8])