wvtf.m

柯恩类时频分布演示例程

function [fa,xieboa,tfr,rtfr,t,f]=wvtf(sig,fs)
N=length(sig);
a=1:1:N
t=20/fs*(a-1)./1000
fc=sig(1:N-100)
fb=sig(1+fs/3:N-100+fs/3)
fa=sig(1+2*fs/3:N-100+2*fs/3)
%subplot(3,1,1);
%plot(fa);
%grid on;
%title('原始信号');
%xlabel('n')
%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,3);
%plot(xieboa);
%title('扰动分量');
%xlabel('n')
%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=xieboa'
sig=(sig-mean(sig))/std(sig,1);
sig=hilbert(sig);
%平滑伪分布
%g=window(@gausswin,M/10+1)
%h=window(@gausswin,M/4+1)
%[tfr,t,f]=tfrspwv(sig,1:M,M,g,h);
M=length(xieboa)
%t=1:M;
%f=0:0.5/M:0.5-0.5/M;
[tfr,t,f]=tfrspwv(sig,1:M,M);
[tfr,rtfr,hat]=tfrrspwv(sig,1:M,M);

%figure(2)
%subplot(211)
%contour(t,f*fs*50,abs(tfr),16);
%title '平滑伪wigner分布'
%xlabel 'n'
%ylabel 'f'
%mesh(t,f,abs(tfr))
%figure(3)
%subplot(212)
%contour(t,f*fs*50,abs(rtfr),16);
%title '重排平滑伪wigner分布'
%xlabel 'n'
%ylabel 'f'
%figure(3)
%mesh(t,f,abs(rtfr));
%title '重排平滑伪wigner分布(三维）'
%xlabel 'n'
%ylabel 'f'