Here is an idea of a partial answer; cf. the answer by Christian Remling. Assume that $(1)$ holds for all $\theta$ in (say) a right neighborhood of $0$ (including $0$), where $|f(x)|=O(1/(1+|x|^k))$ for each $k>0$ (note that the latter condition is not satisfied by the function $f=I_{(0,\infty)}-c$ in the example in the question). Then one can differentiate the integral in $(1)$ in $\theta$ to the right of $0$ to get 
\begin{equation}
	\left. \int_{\mathbb R}\frac{\partial^k g_\theta(x)}{\partial\theta^k} \right|_{\theta=0}\,f(x)\,dx=0,  \tag{2}
\end{equation}
for all $k=0,1,\dots$. 
Expanding $g_\theta(x)/g_0(x)=e^{(\theta x)-(\theta x)^2/2}$ into the Maclaurin series, we have 
\begin{equation*}
	\left. \frac{\partial^k g_\theta(x)}{\partial\theta^k} \right|_{\theta=0}=a_k x^k\,e^{-1/2},   
\end{equation*}
where
\begin{equation*}
	a_k:=\left(-\frac{1}{2}\right)^k k! \sum _{n=\left\lceil k/2\right\rceil }^k
	\frac1{n!}\,\binom{n}{k-n} (-2)^n. 
\end{equation*}
The sequence $(a_k)$ is the sequence A001464 at http://oeis.org/A001464. 

It appears that $a_k\ne0$ for all $k=0,1,\dots$ except $k=2$. If that is indeed so, then, by $(2)$, we would have 
\begin{equation}
	\int_{\mathbb R}x^kf(x)\,dx=0 \tag{3}
\end{equation}
for all $k=0,1,\dots$ except $k=2$. Letting $y:=x^4$, we see that $x^2=\sqrt y$ can be approximated arbitrarily closely by polynomials in $y=x^4$. So, one would have $(3)$ for all $k=0,1,\dots$, including $k=2$, if $|f(x)|$ decreases fast enough as $|x|\to\infty$. This would imply that $f=0$ almost everywhere.