在 Matlab 中对 Wigner-Ville 分布应用窗函数答案

【问题标题】：To apply window function on Wigner-Ville Distribution in Matlab在 Matlab 中对 Wigner-Ville 分布应用窗函数
【发布时间】：2013-12-25 22:25:38
【问题描述】：

我们在考虑 here 如何创建重叠 64 的 Hamming-64 窗口。这是由

h = hamming(64);
h2 = hamming(38);
h = conv(h, h2);

现在，我们正在考虑如何将此窗口函数应用于 Time-Frequency Toolbox 中 Auger 等人的 Wigner-Ville 分布函数的结果变量。函数tfrwv.m没有窗口函数的参数。

所以我们有这些变量

[B,T,F] = tfrwv(data, 1:length(data), length(data));

Here 是相关问题的一个答案，但并不完全相同。有人说将窗口函数应用于结果

逐点相乘

h 的尺寸是 101x1 双倍，而 T 和 F 是 5001x1 双倍。因此，如果逐点相乘，似乎需要对窗口向量进行外推。

再解释here

大约在第二个代码块进行到一半时，我应用了一个窗口缓冲信号的功能。这实际上是一个向量窗口函数与每个缓冲时间块的乘积系列数据。我只是使用偷偷摸摸的对角矩阵技巧来做到这一点高效。

如何将窗口函数应用于变量B、T、和F？

【问题讨论】：

标签： matlab signal-processing time-frequency

【解决方案1】：

查看this paper 关于 Wigner 分发的信息。从第8页到第11页。我认为代码中的tfr(indices,icol) = x(ti+tau,1) .* conj(x(ti-tau,xcol));实现了公式（23）。和和指数部分等价于tfr= fft(tfr);。或者您也可以将我引用的那两行代码视为公式（24）。注意：在论文和代码中，他们假设信号是周期性的N/2，其中N 是length(data)。没关系，您不需要在此处更改您的data。他们只是对原始数据进行了某种扩展。

引自论文，在处理 WDF 之前，将修改后的汉明窗应用于时域信号，以减少由不连续性引起的泄漏数据的有限记录，这将被称为数据逐渐减少。 据我了解，您可以在这里做的是

data1 = conv(h,data);
[B,T,F] = tfrwv(data1, 1:length(data1), length(data1));

我的回答是根据您的实现完成的。你现在可以试一试。

我不完全理解的是您创建重叠 60% 的 Hamming-64 窗口的方法。在spectrogram 中，代码将您的data 拆分为每个长度为64 的小段。如果您想达到spectrogram 与tfrwv 相同的效果，我想您可能还需要拆分您的@987654335 @, 并分别使用conv(data(1:64),hamming(64)), conv(data(38:101),hamming(64)),conv(data(76:139),hamming(64)),....作为tfrwv的输入。

【讨论】：

请看我对你答案的扩展。我们现在知道我最初对 hamming-64 的尝试是不正确的。让我们试试你的解决方案。当我尝试将输入用于 tfrwv 时，我收到错误 dl.dropboxusercontent.com/u/62073194/…。
输入为 data = conv(data(1:64),hamming(64)); [B,T,F] = tfrwv(数据, 1:长度(数据), 长度(数据));然后你用输入数据重新运行程序 = conv(data(38:101),hamming(64)); [B,T,F] = tfrwv(data, 1:length(data), length(data));.... 以此类推。你可以尝试几个片段来观察节拍是否合理，然后编写一个循环对所有片段进行操作
我上一条评论有错字，应该是 data1 = conv(data(1:64),hamming(64)); [B,T,F] = tfrwv(data1, 1:length(data1), length(data1));仍然抛出错误？
我终于找到了我下面代码中的错误。再次感谢您的耐心！我在这里提出了一个关于我当前问题的新问题stackoverflow.com/questions/20860895/…

【解决方案2】：

对 lennon310 答案的扩展。

我跑

data1 = conv(h,data);
[B,T,F] = tfrwv(data1, 1:length(data1), length(data1));

然后用没有窗口的图片来绘制它

我知道没有窗口的图片是正确的，因为我可以看到 60 Hz 水平线是由 MIT-BIH 心律失常数据库中的测量系统设置引起的。患者68岁，有陈旧性心肌梗塞，每秒心跳65次，适合这种方法。

我的原始版本 hamming-64 的实现图片不正确。每分钟跳动这样的人不会活很久。

【讨论】：

【解决方案3】：

lennon310 答案的第二次扩展

data1 = conv(data(1:64),hamming(64)); [B,T,F] = tfrwv(data1, 1:length(data1), length(data1));    
data1 = conv(data(38:101),hamming(64)); [b,t,f] = tfrwv(data1, 1:length(data1), length(data1));

我不知道我应该如何组合从 b 到 B 的数据片段。矩阵对角化似乎不是合适的选择。

为了简化，我只使用hamming(64) 运行，讨论了here 关于矩阵对角化的实现

    B = 0; T = 0; F = 0;    
    data1 = conv(data(1 : 64),hamming(64)); [B,T,F] = tfrwv(data1, 1:length(data1), length(data1));    
    for i=1:10
        data1 = conv(data( 1 + i*37 : 64 + i*37 ),hamming(64)); [b,t,f] = tfrwv(data1, 1:length(data1), length(data1));    
        B = blkdiag(B,b); T = [T t]; F = [F; f];
    end

我明白了

这不是正确的结果。问题在于理解矩阵 B 应该是什么。 将b添加到B后的矩阵B应该是什么样子？

【讨论】：

我认为 size(B(1:floor(numel(F)/2), :))) 应该和 (numel(T),numel(F)) 一致？
在执行B = [B b]之前，还要确认T和t、B和b、F和f的大小一致； T = [T t]； F = [F f]；（1*64，或 64*1 等）
你知道矩阵B包含哪些点吗？我认为这是时间频率的数据。所以想想每次迭代后垂直和水平附加是否合适。
在你的循环中 size(b,1)=127？如果是这样，b 的行就是时间。使用 B=[B b'];此外，使用 filter 代替 conv： filter(data(1 : 64),1,hamming(64))
@lennon310 使用 filter 命令时 b 的大小现在为 64。使用卷积命令时为 127。我在下面添加了您答案的第三个扩展名。

【解决方案4】：

在我对 lennon310 答案的第三次扩展中有一个错误及其症状。 lennon310 答案的第四次扩展

我跑

h = hamming(64);
h2 = hamming(38);
h = conv(h, h2);    
B = 0; T = 0; F = 0;    
data1 = filter(data(1 : 64),1,h); [B,T,F] = tfrwv(data1, 1:length(data1), length(data1));    
for i=1:133
    data1 = filter(data( 1 + i*37 : 64 + i*37 ),1,h); [b,t,f] = tfrwv(data1, 1:length(data1), length(data1));    
    B = [B b']; T = [T t]; F = [F; f];       
end
data1 = filter(data(4959 : 5001),1,h); [b,t,f] = tfrwv(data1, 1:length(data1), length(data1));    
B = [B b']; T = [T t]; F = [F; f];       
T = 49.8899*T;      % dummy constant to get appropriate time interval

得到这样的图片

我没有设法在一张图片中显示所有细峰。一个关于它的新问题here。

我正在策划这个

t = 1/360; % 360 samples per second
fs = 360.5;
imagesc(T*t, F*fs, abs(B))

该算法正在将点累积到正确的维度。 我不确定乘以虚拟常数是否是更早的正确方法。

【讨论】：