I'm certainly not an expert on this subject, but I do have a few naive suggestions.
The function $f(x,y) = \sum_{n,m=-N}^N t_{mn}e^{inx+imy}$ is just a finite Fourier sum defined on a square. It is well-known how to find a finite Fourier sum that "best" approximates a given continuous function on a square.
Therefore, one possible approach would be the following:
Choose a continuous function $F(x,y)$ on the square whose zero set is the given fractal.
Find the Fourier sum $f$ which best approximates the function $F$.
For step (1), one potential choice would be to use standard distance function to a compact set, i.e. $$ F(p) \;=\; \min\{ d(p,q) \mid q\in\text{the fractal set} \} $$ where $d$ denotes Euclidean distance in the plane. The problem with this choice is that, if you perturb the function $F$ slightly, the zero set might disappear entirely.
It would be better to start with a continuous function $F$ whose graph intersects the $xy$-plane transversely along the fractal curve. For something like the the Koch snowflake, you could use a continuous function $F$ which is positive outside of the snowflake and negative inside, e.g. use the distance function for the outside and the negative of the distance function on the inside.
Also, practically speaking, there's no reason why the function $F$ really needs to be continuous. For example, if you start with a piecewise function $F$ which is $1$ outside of the Koch snowflake and $-1$ inside the Koch snowflake, then the Fourier approximations for $F$ might work fairly well.
Edit: I tried this in Mathematica, and it seems to work. Here is a plot of the zero set for a Fourier series with $2500$ terms ($N = 50$):N=50$:
For $F$, I used a function which is $1$ outside the third iterate for the Koch snowflake, and $-1$ inside. I used Green's Theorem to compute the double integrals for the coefficients.
