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} & \mbox{row with 1} \\ 00 & 0(0 \oplus 1) = 01 \\ 01 & 0(1 \oplus 1) = 00 \\ 10 & 1(1 \oplus 1) = 10 \\ 11 & 1(1 \oplus 0) = 11 \end{array} $ This gives us the matrix $ \left( \begin{array}{cc} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{array} \right) $ This generalizes to when $y$ is more than a single bit.