Given two $n$-variable boolean functions $f$ and $g$, we study the problem of computing an $\varepsilon$-approximate isomorphism between them. I.e.\ a permutation $\pi$ of the $n$ variables such that $f(x_1,x_2,\ldots,x_n)$ and $g(x_{\pi(1)},x_{\pi(2)},\ldots,x_{\pi(n)})$ differ on at most an $\varepsilon$ fraction of all boolean inputs $\{0,1\}^n$. We give a randomized $2^{O(\sqrt{n}\log(n)^{O(1)})}$ algorithm that computes a $\frac{1}{2^{\log(n)^{O(1)}}}$-approximate isomorphism between two isomorphic boolean functions $f$ and $g$ that are given by depth $d$ circuits of $n^{O(1)}$ size, where $d$ is a constant independent of $n$. In contrast, the best known algorithm for computing an exact isomorphism between $n$-ary boolean functions has running time $2^{O(n)}$ \cite{Luks99} even for functions computed by $n^{O(1)}$ size DNF formulas. Our algorithm is based on a recent result for hypergraph isomorphism with bounded edge size \cite{BabaiC08} and the classical Linial Mansour-Nisan result on approximating small depth and size boolean circuits by small degree polynomials using Fourier analysis.