平滑伪Wigner-Ville分布计算程序

mengfuse

% Wigner-Ville分布计算程序
function[WVD,t,f]=SWWignerVille(Sig,SampFreq,FreqBins,GLen,HLen)
%  Sig      : 输入信号
%  FreqBins : 频率轴划分区间数（默认值：512）
%  SampFreq : 信号的采集频率
%  alfa     : 窗函数控制参数
%  GLen     : 窗函数g的长度（默认值：FreqBins/5）
%  HLen     : 窗函数h的长度（默认值：FreqBins/4）
%  WVD      :   SWWigner-Ville计算结果

if (nargin<1),
    error('At least 1 input required');
end;
% Sig=hilbert(real(Sig));    % 计算信号的解析信号
SigLen=length(Sig);        % 获取信号的长度
if (nargin<2),
    SampFreq=1;
elseif (nargin<3),
    FreqBins=512;
elseif (nargin<4),
    GLen=floor(FreqBins/5);   
    HLen=floor(FreqBins/4);
end;
FreqBins=2^nextpow2(FreqBins);
FreqBins=min(FreqBins,SigLen);
GLen=GLen+1-rem(GLen,2);
HLen=HLen+1-rem(HLen,2);
GWin=window(@gausswin,GLen);
HWin=window(@gausswin,HLen);
Lg=(GLen-1)/2;
Lh=(HLen-1)/2;
HWin=HWin/HWin(Lh+1);
WVD=zeros(FreqBins,SigLen);
wait=waitbar(0,'Under calculation,please wait...');

for kk=1:SigLen,
    waitbar(kk/SigLen,wait);
    MTau=min([kk+Lg-1,SigLen-kk+Lg,round(FreqBins/2)-1,Lh]);
    k=-min([Lg,SigLen-kk]):min([Lg,kk-1]);
    SubG=GWin(Lg+1+k);
    SubG=SubG/sum(SubG);
    WVD(1,kk)=sum(SubG.*Sig(kk-k,1).*conj(Sig(kk-k)));  %
    for tau=1:MTau,
      k=-min([Lg,SigLen-kk-tau]):min([Lg,kk-tau-1]);
      SubG=GWin(Lg+1+k);
      SubG=SubG/sum(SubG);
      R=sum(SubG.*Sig(kk+tau-k,1).*conj(Sig(kk-tau-k)));
      WVD(1+tau,kk)=HWin(Lh+tau+1)*R;
      R=sum(SubG.*Sig(kk-tau-k,1).*conj(Sig(kk+tau-k)));
      WVD(FreqBins+1-tau,kk)=HWin(Lh-tau+1)*R;
    end;   
end;
close(wait);
WVD=fft(WVD);
f=linspace(0,0.5,FreqBins)*SampFreq;
t=(0:(SigLen-1))/SampFreq;
set(gcf,'position',[20 100 500 430]);
set(gcf,'color','w');
axes('position',[0.1 0.45 0.53 0.5]);
mesh(t,f,abs(WVD));
axis([min(t) max(t) min(f) max(f)]);
colorbar;
xlabel('t/s');
ylabel('f/Hz');
title('平滑伪Wigner-Ville分布');

输入以下：
%  SampFreq=100;
%  FreqBins=512;
%  t=0:1/SampFreq:4;
%  Sig=sin(20*pi*t)+sin(80*pi*t);
%  WignerVille(Sig,SampFreq,FreqBins,Freqbins/5,FreqBins/4);

出现以下错误：
??? Index exceeds matrix dimensions.

Error in ==> SWWignerVille at 43
    WVD(Temp,kk)=sum(SubG.*Sig(kk-k,1).*conj(Sig(kk-k)));

谁能给修改下。谢谢

mengfuse

犟牛

1. WVD(1,kk)=sum(SubG.*Sig(kk-k,1).*conj(Sig(kk-k)));

复制代码

改为：

1. WVD(1,kk)=sum(SubG'.*Sig(kk-k).*conj(Sig(kk-k)));

复制代码

后面还有类似错误，自己改吧

ChaChing

Ref: 12F, 常见的程序出错问题整理 http://forum.vibunion.com/thread-46001-1-1.html
from http://forum.vibunion.com/home-s ... -blog-id-18250.html

tf19830910

想请教楼主一下，这个平滑算法是根据什么文献编写的？时频工具箱中的算法似乎和你是一样的，但是我找不到这个算法的理论文献。