Multi-qubit gates matrix representation
The first thing to recognize is that this is a permutation matrix where each row and column contains a single $1$ entry with the remaining entries being $0$.  This follows from the fact that $x$ and $y$ represent classical bit strings and $f$ is a classical boolean function.  

Let's take the basic case where $x$ is from the space of $n$ bit strings, $y$ is a single bit, and $\mathcal{f}$ maps $\{0,1\}^n$ to $\{0,1\}$. The operation then functions on $n+1$ bits and is represented by a $2^{n+1} \times 2^{n+1}$ matrix. If you let the binary number $xy$ index the matrix columns, then for column $xy$ set row $x(f(x) \oplus y)$ to 1 and all other rows to 0. 

Let's look at an example from Deutsch's algorithm. Let $x$ be from $\{0,1\}$ and $f$ be bit negation, i.e. $f(0) = 1$ and $f(1) = 0$.  From here we compute, for each row, the column number containing the 1.

$\begin{array}{cc}
\mbox{column} &amp; \mbox{row with 1} \\
00 &amp; 0(0 \oplus 1) = 01 \\
01 &amp; 0(1 \oplus 1) = 00 \\
10 &amp; 1(1 \oplus 1) = 10 \\
11 &amp; 1(1 \oplus 0) = 11 
\end{array} $

This gives us the matrix

$ \left( \begin{array}{cc}
0 &amp; 1 &amp; 0 &amp; 0 \\
1 &amp; 0 &amp; 0 &amp; 0 \\
0 &amp; 0 &amp; 1 &amp; 0 \\
0 &amp; 0 &amp; 0 &amp; 1 
\end{array} \right) $

This generalizes to when $y$ is more than a single bit.