choi_willaims.m

滤波器组的PWVD很不错的我已经试过了很好

function [wx,t_lable,f_lable] = Choi_Willaims (x0,nfft,sigma)
%计算 Choi-Williams 时频分布
%	x0     - 信号序列
%	nfft  - 进行FFT计算的长度。要求比信号的2倍长度长，且为2 的幂次       
%	sigma - sigma for the Choi-Williams window  [default = 0.05]
% wx :Choi-Williams 分布
%t_lable,f_lable :为时频图的时间轴和频率轴

lx = length(x0);
lfft = 2^nextpow2(2*lx);                      % minimum FFT length
if (exist('nfft') ~= 1) nfft = lfft; end
if (isempty(nfft)) nfft = lfft; end

if (nfft < 2*lx)
   disp(['FFT length must exceed twice the signal length'])
   disp(['     resetting FFT length to ',int2str(lfft)])
   nfft = lfft;
end

x  = zeros(nfft,1);  x(1:lx) = x0(:);       cx = conj(x);
wx = zeros(nfft,nfft);
L1 = lx - 1;


for n=0:L1
   indm = max(-n,-L1+n) : min(n,L1-n);
   indy = indm + (indm < 0) * nfft ;   
   y = zeros(nfft,1);
   y(indy + 1) = x(n+indm + 1) .* cx(n-indm + 1);
   wx(:,n+1)   = y;
end


wx  =fftshift ( ifft(wx.').' );
if (finite(sigma))
     win =   [ (1:nfft)-nfft/2-1 ]' * [ (1:nfft)-nfft/2-1] ...
     		/ nfft;
     win = ( exp (- win.^2 / sigma) );
else
     win = ones(nfft,nfft);
end
wx  = wx .* win;
wx  = fft2(wx); % fft along both time and frequency
wx = abs(wx) ;


wx  = wx(:,1:length(x0));
[m ,n]= size(wx);
t_lable = 1:n ; f_lable = (1:m)/(2*nfft) ;
contour(t_lable,f_lable ,wx,4) ;
xlabel('time ') ;ylabel('frequency')
title(['choi willians distributrion ' ' sigma=',num2str(sigma)]);
zoom on;