Let me change your notation just a little to make answering easier. $y_{im}$ is the value of the outcome for experiment $i$ at frequency $m$, $i=1,\ldots,N$ and $m=1,\ldots,M$. You have a function, $h(x,\theta)$, which outputs a $M$-dimensional vector of predicted $y$, one for each frequency, given a particular value for predictor variables, $x$, and parameter value $\theta$. Connecting to your notation, $h(x,\theta)=g(f,x,\theta)$ (I think). It's going to be convenient to think of a whole experiment's observations at once, so let $Y_i$ be the $M$-vector of all of experiment $i$'s $y$ for the various frequencies stacked up.
The data-generating process is then hypothesized to be:
\begin{align}
Y_i &= h(x_i,\theta)+\epsilon_i
\end{align}
Measurement error (or whatever the source of error is) for experiment $i$ is described by the random variables $\epsilon_i$. All of the objects in the above equation are $M$-vectors. You say that different frequencies' data are likely to be correlated. I assume it is also the case that different frequencies' data could have different variances. If we are safe to assume that these variances and covariances are the same across experiments (and that each experiment is independent of each other experiment), then we can assume that $V(\epsilon_i)=\Sigma$, where $\Sigma$ is an $M \times M$ variance-covariance matrix, and $V(\epsilon)=I \otimes \Sigma$, where $\epsilon$ is the "stacked-up" vector of errors ($NM$ by 1) and $\otimes$ is Kronecker product.
The above set-up is called non-linear seemingly unrelated regression. The standard reference in Economics is Gallant, AR (1975) Seemingly unrelated nonlinear regressions. Journal of Econometrics, 3(1): 35-50.
The way you estimate the parameters (if you know $\Sigma$) is by solving:
\begin{align}
min_\theta \sum_{i=1}^N \left( Y_i-h(X_i,\theta)\right)' \Sigma^{-1} \left( Y_i-h(X_i,\theta)\right)
\end{align}
The variance of the (consistent and asymptotically normal) estimator, $\hat{\theta}$, so-defined is:
\begin{align}
V(\hat{\theta}) &= \left( \sum_{i=1}^N \frac{\partial h}{\partial \theta}' \Sigma \frac{\partial h}{\partial \theta} \right)^{-1}
\end{align}
Normally, you don't know what $\Sigma$ is, however, so you have to estimate it. The way you do this is by picking a trial value of $\Sigma$ (like the identity matrix). Then you do the minimization above to get a first try at an estimator of $\theta$, call it $\tilde{\theta}$. Then you "estimate" the $\epsilon$ by the residuals, $e_i=Y_i-h(X_i,\tilde{\theta})$ and calculate their sample variance matrix:
\begin{align}
\hat{\Sigma} &= \frac{1}{N} \sum_{i=1}^N e_i e_i'
\end{align}
Finally, you use $\hat{\Sigma}$ in place of $\Sigma$ above in the minimization and the variance-covariance matrix. This procedure is called two-stage, non-linear seemingly unrelated regression (or words to that effect).
In R, I think the function you want is called nlsystemfit which has documentation here.
