dqfenjie.m

对信号虚构三相

clear all;
hold off;
N=800;
fs=100;
x=1:1:N;%总共取4个周期的波形，800点
for a=1:1:length(x)
    t(a)=20/fs*(x(a)-1)/1000
   if a<300
       f(a)=sin(2*50*pi*t(a))
   else if a>600
           f(a)=sin(2*50*pi*t(a));
       else
       f(a)=0.0*sin(2*50*pi*t(a))+0.0*sin(6*50*pi*t(a))+0.0*sin(10*50*pi*t(a));
       
       end
   end
end
fm1=fmconst(800,0.3);
am1=amexpo1s(800,300,40)
r1=1*am1.*fm1;
for a=1:1:length(x)
    if t(a)<0.05
        r(a)=0;
    else if t(a)>0.07
            r(a)=0;
        else
            r(a)=-0.0*sin(2*50*pi*t(a))
        end
    end
end
f=f+0.0*randn(1,length(f))+0*real(r1')
fc=f(1:600);
fb=f(1+fs/3:600+fs/3)
fa=f(1+2*fs/3:600+2*fs/3)
subplot(3,1,1);
plot(fc);
grid on;
title('（a）电压暂降信号');
xlabel('采样点')
ylabel('U/kv');
subplot(3,1,2);
hold on;
%plot(fb);
grid on;
%title('b相');
hold on;
%subplot(3,1,3);
%plot(fc);
%grid on;
%title('c相');
xlabel 't/s';
ylabel 'U/kV';
hold on;
N=length(fa)
%构造ab分量
C32=(2/3)^0.5*[1 -1/2 -1/2;0 3^0.5/2 -3^0.5/2];
A=[fa;fb;fc];
X=C32*A;

%gouzaodq分量
for n=1:1:N
ua(n)=X(1,n);
ub(n)=X(2,n);
  B=[sin(2*pi*50*t(n)) -cos(2*pi*50*t(n));-cos(2*pi*50*t(n)) -sin(2*pi*50*t(n))];
  C=[ua(n);ub(n)]
  X1=B*C;
    ud(n)=X1(1,1);
    uq(n)=X1(2,1);
    U2(n)=ud(n)^2+uq(n)^2;
end;
%quzao1
M=60;
lg=2;
input=ud;
ud1=quzao1(ud,60,2,64,2);
uq1=quzao1(uq,60,2,64,2);
t1=t(1:length(ud1));
%hold off;
%plot(t1,ud1,'r');
N=length(ud1);
ud2=quzao1(ud1,60,2,64,2);
uq2=quzao1(uq1,60,2,64,2);
%hold off;
%plot(t2,uq2,'r');
%dq分量fanbianhuan
N=length(ud2);
for n=1:1:N
  B=[sin(2*pi*50*t(n)) -cos(2*pi*50*t(n));-cos(2*pi*50*t(n)) -sin(2*pi*50*t(n))];
  Bn=inv(B);
  X2=Bn*[ud2(n);uq2(n)];
    ua1(n)=X2(1,1);
    ub1(n)=X2(2,1);
end;
%a bfanbianhaun
for n=1:1:N
 D=[ua1(n);ub1(n)];
 C23=(2/3)^0.5*[1 0;-1/2 3^0.5/2;-1/2 -3^0.5/2];
   X3=C23*D;
   faj(n)=X3(1,1);
   fbj(n)=X3(2,1);
   fcj(n)=X3(3,1);
end;
%hold off;
%subplot(3,1,2);
%plot(faj);
%grid on;
%title('基波分量');
%xlabel('n')
%ylabel('U/kv');
grid on;
hold on;
%subplot(3,1,2);
%plot(t,fbj,'r');
hold on;
%subplot(3,1,3);
%plot(t,fcj,'k');
hold on;
xieboa=fa-faj;
xiebob=fb-fbj;
xieboc=fc-fcj;
subplot(3,1,2);
plot(xieboc);
title('（b）扰动分量');
xlabel('采样点')
ylabel('U/kv');
%grid on;
hold on;
%subplot(3,1,2);
%plot(t,xiebob,'r');
hold on;
%subplot(3,1,3);
%plot(t,xieboc,'k');
hold on;
%%%%%%%%%%%%%%%%%%%%
%信号频率归一化和hilbert变换
sig=xieboc'
sig=(sig-mean(sig))/std(sig,1);
sig=hilbert(sig);
%平滑伪分布
g=window(@Hamming,51)
%h=window(@Hamming,151)
%[tfr,t,f]=tfrspwv(sig,1:M,M,g,h);
M=length(xieboa)
[tfr,t,f]=tfrspwv(sig,1:M,M,g);
[tfr,rtfr,hat]=tfrrspwv(sig,1:M,M,g);
%[tfr,t,f]=tfrspwv(sig,1:M,M);
%figure(2)
%subplot(211)
%contour(t,f*5000,abs(tfr),16);
%title '（a）平滑伪wigner分布'
%xlabel '采样点'
%ylabel 'f/Hz'
%mesh(t,f,abs(tfr))
%figure(3)
subplot(313)
contour(t,f*5000,abs(rtfr),16);
%mesh(t,f*5000,abs(rtfr));
title '（c)重排平滑伪wigner-vill分布'
xlabel '采样点'
ylabel 'f/Hz'
%计算边缘特性
for i=1:1:600
    yy1(i)=sum(tfr(i,1:600));
end
%计算边缘特性
for i=1:1:600
    yy2(i)=sum(rtfr(i,1:600));
end
figure(2);
%subplot(121)
%plot(f*5000,yy1);
%subplot(122)
plot(f*5000,yy2);
ylabel ‘能量强度’
%xlabel ‘f/Hz'