二维矩阵中多个一维信号与二维矩阵中多个一维核的卷积

Question

我有一个随机定义的 H 大小为 600 x 10 的矩阵。此矩阵 H 中的每个元素都可以表示为 H(k,t)。我获得了一个语音频谱图 S，它是 600 x 597。我使用 Mel 功能获得它，所以它应该是 40 x 611 但后来我使用了帧堆叠概念，其中我将 15 帧堆叠在一起。因此它给了我 (40x15) x (611-15+1) 即 600 x 597。

现在我想获得一个输出矩阵Y，它由基于卷积的方程Y(k,t) = ∑ H(k,τ)S(k,t-τ)给出。总和从 τ=0 到 τ=Lh-1。 Lh 在这种情况下为 597。

不知道如何获取Y。另外，我怀疑在计算卷积时对 H 和 S 的索引。具体来说，对于 Y(1,1)，我们有：

Y(1,1) = H(1,0)S(1,1) + H(1,1)S(1,0) + H(1,2)S(1,-1) + H(1,3)S(1,-2) + ...

现在，MATLAB 中没有负索引这样的东西 - 例如，S(1,-1) S(1,-2) 等等。那么，我应该使用什么类型的卷积来获得Y呢？我尝试使用 conv2 或 fftfilt，但我认为这不会给我 Y，因为 Y 也必须是 S 的大小。

Answer 1

这很容易。这是仅应用于一维的 2D 信号的卷积。如果我们假设变量 k 用于访问行而 t 用于访问列，则可以将 H 和 S 的每一行视为单独的信号S 的每一行是一维信号，H 的每一行是卷积核。

有两种方法可以解决这个问题。

时域

如果你想坚持使用时域，最简单的方法是遍历输出的每一行，找到 S 和 H 的每一对行的卷积并存储相应输出行中的输出。据我所知，没有任何实用程序可以仅在给定 N-D 信号的情况下在一维上进行卷积……除非您进入频域，但让我们稍后再讨论。

类似于：

Y = zeros(size(S));
for idx = 1 : size(Y,1)
    Y(idx,:) = conv(S(idx,:), H(idx,:), 'same');
end

对于输出的每一行，我们用一行 S 和一行 H 执行逐行卷积。我使用 'same' 标志，因为输出应该与一行 S 的大小相同......这是更大的行。

频域

您也可以在频域中执行相同的计算。如果您对卷积和傅里叶变换的性质有所了解，您就会知道时域中的卷积是频域中的乘法。您对两个信号进行傅里叶变换，将它们按元素相乘，然后再进行傅里叶逆变换。

但是，您需要牢记以下复杂问题：

进行全卷积意味着输出信号的最终长度为length(A)+length(B)-1，假设A和B是一维信号。因此，您需要确保 A 和 B 都是 zero-padded 以便它们匹配相同的大小。 为什么 确保信号大小相同的原因是允许乘法运算。
一旦你在频域中乘以信号然后取反，你会看到 Y 的每一行都是 full 长度的卷积。为确保获得与输入大小相同的输出，您需要 trim 关闭 开头和结尾的一些点。具体来说，由于 H 的每个内核/列长度为 10，因此您必须删除输出中每个信号的前 5 个点和后 5 个点，以匹配您在 for 循环代码中获得的内容。
由于FFT算法的特性，通常在傅里叶逆变换后，会残留一些复系数。最好使用 real 删除结果的复数值部分。

将所有这些理论放在一起，这就是代码的样子：

%// Define zero-padded H and S matrices
%// Rows are the same, but columns must be padded to match point #1
H2 = zeros(size(H,1), size(H,2)+size(S,2)-1);
S2 = zeros(size(S,1), size(H,2)+size(S,2)-1);

%// Place H and S at the beginning and leave the rest of the columns zero
H2(:,1:size(H,2)) = H;
S2(:,1:size(S,2)) = S;

%// Perform Fourier Transform on each row separately of padded matrices
Hfft = fft(H2, [], 2);
Sfft = fft(S2, [], 2);

%// Perform convolution
Yfft = Hfft .* Sfft;

%// Take inverse Fourier Transform and convert to real
Y2 = real(ifft(Yfft, [], 2));

%// Trim off unnecessary values
Y2 = Y2(:,size(H,2)/2 + 1 : end - size(H,2)/2 + 1);

Y2应该是卷积后的结果，应该匹配前面for循环代码中的Y。

两者的比较

如果你真的想比较它们，我们可以。我们首先需要做的是定义 H 和 S。为了重建我所做的，我使用已知种子生成了随机值：

rng(123);
H = rand(600,10);
S = rand(600,597);

一旦我们运行上面的时域版本和频域版本的代码，让我们看看它们在命令提示符中是如何匹配的。让我们显示前 5 行和 5 列：

>> format long g;
>> Y(1:5,1:5)

ans =

          1.63740867892464          1.94924208172753          2.38365646354643          2.05455605619097          2.21772526557861
          2.04478411247085          2.15915645246324          2.13672842742653          2.07661341840867          2.61567534623066
         0.987777477630861           1.3969752201781          2.46239452105228          3.07699790208937          3.04588738611503
          1.36555260994797          1.48506871890027          1.69896157726456          1.82433906982894          1.62526864072424
          1.52085236885395          2.53506897420001          2.36780282057747          2.22335617436888          3.04025523335182

>> Y2(1:5,1:5)

ans =

      1.63740867892464          1.94924208172753          2.38365646354643          2.05455605619097          2.21772526557861
      2.04478411247085          2.15915645246324          2.13672842742653          2.07661341840867          2.61567534623066
     0.987777477630861           1.3969752201781          2.46239452105228          3.07699790208937          3.04588738611503
      1.36555260994797          1.48506871890027          1.69896157726456          1.82433906982894          1.62526864072424
      1.52085236885395          2.53506897420001          2.36780282057747          2.22335617436888          3.04025523335182

我觉得不错！作为另一种衡量标准，让我们找出 Y 中的一个值与 Y2 中的对应值之间的最大差异：

>> max(abs(Y(:) - Y2(:)))

ans =

      5.32907051820075e-15

也就是说，两个输出之间的最大误差约为 10^-15。我会说那很好。

二维矩阵中多个一维信号与二维矩阵中多个一维核的卷积

Convolution of multiple 1D signals in a 2D matrix with multiple 1D kernels in a 2D matrix

matlab

signal-processing

fft

convolution

时域

频域

两者的比较