The reason why you go from 16 to 3 channels is that by convention filters span the entire depth of the input. Therefore, your filters would actually be
16 in order to cover all channels of the input.
The output of the convolution automatically has a depth equal to the number of filters (so in your case this is 3) because you have an (
m x k) filter matrix, where
m is the number of filters and
k is the number of elements in the unrolled filter (in your case,
m = 3 and
k = 7 x 7 x 16 = 784, so the filter matrix is
3 x 784).
The input is usually unrolled according to the im2col procedure, where each tile corresponding to a single filter location is stretched into a column equal to the unrolled filter size. This is repeated for each filter location, so you end up with very large matrix of size (
k x n), where
k is the same as
k above in the filter matrix, and
n depends on your padding and stride.
Multiplying the (
m x k) filter matrix with the (
k x n) input matrix gives you an (
m x n) output matrix, where
m is the number of filters.
You can find some very nice visual explanations of the convolution procedure here and here.