Here's the general derivation of the commonly used, and often (but not always) valid, uncertainty propagation formula for independent small Gaussian errors.
$\newcommand{\bbv}[1]{\mathbf{#1}}$
Consider a quantity $y$, calculated from measured quantities $\bbv{x}$
$$
y + \Delta{y} = f(\bbv{x}+\Delta\bbv{x}) = f(x_1+\Delta{x}_1,x_2+\Delta{x}_2,\ldots,x_n+\Delta{x}_n)
$$
Note that I don't write $\pm$ symbols: it's the actual, positive or negative, deviation of a single measured value from the true values that I mean. The $\Delta x$ are not the measurement uncertainties, but randomly-distributed physical quantities that have somehow found their way into our measurement.
Assume now the deviations to be
- independent, i.e. stemming from measurements not directly coupled
- small
where by small I mean, we can Taylor expand the function:
$$
y + \Delta{y} \approx f(\bbv{x}) + \Delta\bbv{x}\cdot\bbv\nabla{f(\bbv{x})}
= f(\bbv{x}) + \Delta{x_1}\cdot\frac{\partial}{\partial{x_1}}{f} + \Delta{x_2}\cdot\frac{\partial}{\partial{x_2}}{f} + \ldots
$$
with only the terms I've writte there. Since $f(\bbv{x})=y$, it's just
$$
\Delta{y} = \Delta{x_1}\cdot\partial_1{f} + \Delta{x_2}\cdot\partial_2{f} + \ldots
$$
(that's not the result yet!)
That the errors are independent means that there exists a probability density function for each of the measured quantities $\psi_i(\Delta{x_i})$. What we're interested in is the corresponding PDF of the result uncertainties
$\newcommand{\ttd}{\mathrm{d}}$
$$\begin{aligned}
\phi(\Delta{y}) =& \text{probability density of }\bigl(f(\bbv{x}+\Delta\bbv{x})=y+\Delta{y}\bigr)
\\=& \int_{\{(\Delta x_1,\ldots,\Delta x_n)\in\mathbb{R}^n:f(\bbv{x}+\Delta\bbv{x})=y+\Delta{y}\}}\!\!\!\!\!\!\!\!\!\!\!\!\ttd\!\Delta\!x_1\!\cdot\!\psi_1(\Delta x_1)\ \ttd\!\Delta\!x_2\!\cdot\!\psi_2(\Delta x_2)\ \ldots
\end{aligned}$$
integrating over an $(n-1)$-dimensional manifold (this strictly needs to be written in a less ambiguous way, but it's ok for our purpose).
Now put in the Taylor series approximation
$$\begin{aligned}
\phi(\Delta{y}) =& \int_{\{(\Delta x_1,\ldots,\Delta x_n):\Delta{y} = \Delta{x_1}\cdot\partial_1{f} + \Delta{x_2}\cdot\partial_2{f}+\ldots\}}\!\!\!\!\!\!\!\!\!\!\!\!\ttd\!\Delta\!x_1\!\cdot\!\psi_1(\Delta x_1)\ \ttd\!\Delta\!x_2\!\cdot\!\psi_2(\Delta x_2)\ \ldots
\end{aligned}$$
and note that this is just the convolution
$$
(f_1\star f_2\star\ldots\star f_n)(y) = \int_{\{(x_1,x_2,\ldots,x_n)\colon\sum_{x_n}=y\}}\!\!\!\!\!\!\!\!\!\!\!\!\ttd x_1\!\cdot\!f_1(x_1)\ \ttd x_2\!\cdot\!f_2(x_2)\ \ldots
$$
of the functions $\psi_i$. If these are well-behaved in a very weak sense, we can apply the convolution theorem
$$
\bigl(\mathrm{FT}(f_1\star f_2\star\ldots\star f_n)\bigr)(k) = \prod_i\mathrm{FT}(f_i)(k)
$$
(depending on your definition of the Fourier transform, with some scaling factor). In our case,
$$
\mathrm{FT}(\phi)(k) = \prod_i\mathrm{FT}\left(\lambda x.\psi_i\Bigl(\frac{x}{\partial_if}\Bigr)\right)(k)
$$
To actually compute this now, we need to assume some specific shape of the functions $\psi_i$, i.e. of the "shape" of the deviations in our measurements. And for such things, a normal distribution is very often a good bet:
$$
\psi_i(\Delta{x}) = N_i\cdot \exp\Bigl(\frac{-\Delta\!x^2}{2\cdot\sigma{x}_i{}^2}\Bigr).
$$
(Why did nobody correct me? I'd written complete rubbish here!)
We're not really interested in the normalization constant $N_i$, but the standard deviation $\sigma{x}_i$ is highly relevant: this is the actual uncertainty of the measured quantity $x_i$, the thing that's normally meant when writing $\Delta{x_i}$.
The Fourier transform of such a Gaussian bell curve is, conveniently, again a normal distribution:
$$
\mathrm{FT}(\psi_i)(k) \propto \exp\Bigl(\frac{-k^2\cdot \sigma{x}_i{}^2}{2}\Bigr).
$$
We actually need "stretched" versions of the functions,
$$
\mathrm{FT}\bigl(\lambda x.\psi_i(\tfrac{x}{s})\bigr)(k)
= \mathrm{FT}\left(\exp\Bigl(\frac{-\frac{\Delta\!x^2}{s^2}}{2\cdot\sigma{x}_i{}^2}\Bigr)\right)(k)
\propto e^{\frac{-k^2\cdot s^2\cdot\sigma{x}_i{}^2}{2}}.
$$
Using this, we obtain
$$\begin{aligned}
\mathrm{FT}(\phi)(k) \propto \prod_i\exp\Bigl(\frac{-k^2\cdot (\partial_if)^2\cdot \sigma{x}_i{}^2}{2}\Bigr)
\\=& \exp\Bigl(\frac{-k^2}2\sum_i(\partial_if\cdot \sigma{x}_i)^2\Bigr)
\end{aligned}$$
which is yet again a Gaussian function, so we can easily back-transform it,
$$\begin{aligned}
\phi(\Delta{y}) \propto \mathrm{IFT}\left(\exp\Bigl(\frac{-k^2}2\sum_i(\partial_if\cdot \sigma{x}_i)^2\Bigr)\right)(\Delta{y})
\\\propto& \exp\Bigl(\frac{-\Delta\!y^2}{2\sum_i(\partial_if\cdot \sigma{x}_i)^2}\Bigr).
\end{aligned}$$
Here, we can associate $\sum_i(\partial_if\cdot \sigma{x}_i)^2$ with the (squared) standard deviation of this normal distribution. So, summarizing, the uncertainty of the quantity
$$
y = f(\bbv{x})
$$
is
$$
\sigma{y} = \sqrt{\sum_i\Bigl(\frac{\partial f}{\partial x_i}\cdot \sigma{x_i}\Bigr)^2}
$$
and that's it!