× {{alert.msg}} Never ask again
Get notified about new tutorials RECEIVE NEW TUTORIALS

Understanding behaviour of MATLAB's convn

Ray Phan
Feb 02, 2015
<p>You <strong>almost</strong> have it correct. There are two things slightly wrong with your understanding:</p> <ol> <li><p>You chose <code>valid</code> as the convolution flag. This means that the output returned from the convolution has its size so that when you are using the kernel to sweep over the matrix, it has to fit comfortably inside the matrix itself. Therefore, the first "valid" output that is returned is actually for the computation at location <code>(2,2,1)</code> of your matrix. This means that you can fit your kernel comfortably at this location, and this corresponds to position <code>(1,1)</code> of the output. To demonstrate, this is what <code>a</code> and <code>ker</code> look like for me using your above code:</p> <pre><code>&gt;&gt; a a(:,:,1) = 0.9930 0.2325 0.0059 0.2932 0.1270 0.8717 0.3560 0.2365 0.3006 0.3657 0.6321 0.7772 0.7102 0.9298 0.3743 0.6344 0.5339 0.0262 0.0459 0.9585 0.1488 0.2140 0.2812 0.1620 0.8876 0.7110 0.4298 0.9400 0.1054 0.3623 0.5974 0.0161 0.9710 0.8729 0.8327 a(:,:,2) = 0.8461 0.0077 0.5400 0.2982 0.9483 0.9275 0.8572 0.1239 0.0848 0.5681 0.4186 0.5560 0.1984 0.0266 0.5965 0.2255 0.2255 0.4531 0.5006 0.0521 0.9201 0.0164 0.8751 0.5721 0.9324 0.0035 0.4068 0.6809 0.7212 0.3636 0.6610 0.5875 0.4809 0.3724 0.9042 &gt;&gt; ker ker(:,:,1) = 0.5395 0.4849 0.0970 0.3418 0.6263 0.9883 0.4619 0.7989 0.0055 0.3752 0.9630 0.7988 ker(:,:,2) = 0.2082 0.4105 0.6508 0.2669 0.4434 0.1910 0.8655 0.5021 0.7156 0.9675 0.0252 0.0674 </code></pre> <p>As you can see, at position <code>(2,2,1)</code> in the matrix <code>a</code>, <code>ker</code> can fit comfortably inside the matrix and if you recall from convolution, it is simply a sum of element-by-element products between the kernel and the subset of the matrix at position <code>(2,2,1)</code> that is the same size as your kernel (actually, you need to do something else to the kernel which I will reserve for my next point - see below). Therefore, the coefficient that you are calculating is actually the output at <code>(2,2,1)</code>, not at <code>(1,1,1)</code>. From the gist of it though, you already know this, but I wanted to put that out there in case you didn't know.</p></li> <li><p>You are forgetting that for N-D convolution, you need to <strong>flip the mask in each dimension</strong>. If you remember from 1D convolution, the mask must be flipped in the horizontally. What I mean by flipped is that you simply place the elements in reverse order. An array of <code>[1 2 3 4]</code> for example would become <code>[4 3 2 1]</code>. In 2D convolution, you must flip both horizontally and vertically. Therefore, you would take each row of your matrix and place each row in reverse order, much like the 1D case. Here, you would treat each row as a 1D signal and do the flipping. Once you accomplish this, you would take this flipped result, and treat each <strong>column</strong> as a 1D signal and do the flipping again.</p> <p>Now, in your case for 3D, you must flip horizontally, vertically <strong>and temporally</strong>. This means that you would need to perform the 2D flipping for each slice of your matrix independently, you would then grab single columns in a 3D fashion and treat those as 1D signals. In MATLAB syntax, you would get <code>ker(1,1,:)</code>, treat this as a 1D signal, then flip. You would repeat this for <code>ker(1,2,:)</code>, <code>ker(1,3,:)</code> etc. until you are finished with the first slice. Bear in mind that we don't go to the second slice or any of the other slices and repeat what we just did. Because you are taking a 3D section of your matrix, you are inherently operating over all of the slices for each 3D column you extract. Therefore, only look at the first slice of your matrix, and so you need to do this to your kernel before computing the convolution:</p> <pre><code>ker_flipped = flipdim(flipdim(flipdim(ker, 1), 2), 3); </code></pre> <p><a href="http://www.mathworks.com/help/matlab/ref/flipdim.html" rel="nofollow"><code>flipdim</code></a> performs the flipping on a specified axis. In our case, we are doing it vertically, then taking the result and doing it horizontally, and then again doing it temporally. You would then use <code>ker_flipped</code> in your summation instead. Take note that it doesn't matter which order you do the flipping. <code>flipdim</code> operates on each dimension independently, and so as long as you remember to flip all dimensions, the output will be the same.</p></li> </ol> <hr> <p>To demonstrate, here's what the output looks like with <code>convn</code>:</p> <pre><code>c = 4.1837 4.1843 5.1187 6.1535 4.5262 5.3253 5.5181 5.8375 5.1311 4.7648 5.3608 7.1241 </code></pre> <p>Now, to determine what <code>c(1,1)</code> is by hand, you would need to do your calculation on the <strong>flipped</strong> kernel:</p> <pre><code>ker_flipped = flipdim(flipdim(flipdim(ker, 1), 2), 3); c11 = sum(sum(a(1:3,1:4,1).*ker_flipped(:,:,1)))+sum(sum(a(1:3,1:4,2).*ker_flipped(:,:,2))); </code></pre> <p>The output of what we get is:</p> <pre><code>c11 = 4.1837 </code></pre> <p>As you can see, this verifies what we get by hand with the calculation done in MATLAB using <code>convn</code>. If you want to compare more digits of precision, use <code>format long</code> and compare them both:</p> <pre><code>&gt;&gt; format long; &gt;&gt; disp(c11) 4.183698205668000 &gt;&gt; disp(c(1,1)) 4.183698205668001 </code></pre> <p>As you can see, all of the digits are the same, except for the last one. That is attributed to numerical round-off. To be absolutely sure:</p> <pre><code>&gt;&gt; disp(abs(c11 - c(1,1))); 8.881784197001252e-16 </code></pre> <p>... I think a difference of an order or 10<sup>-16</sup> is good enough for me to show that they're equal, right?</p> <p>This tip was originally posted on <a href="http://stackoverflow.com/questions/27852571/Understanding%20behaviour%20of%20MATLAB's%20convn/27854391">Stack Overflow</a>.</p>
comments powered by Disqus