Fit isn't restricted to polynomials. You should choose your basis functions to match the nature of the problem. For example, suppose I believe my data has a jump (of unknown height) at time 0, added to an unknown linear trend with an unknown offset and measurement noise. Something like this:

data = Table[{x, 
   3 UnitStep[x] - 0.3 x + 2 + 
    RandomVariate[NormalDistribution[]]}, {x, -10, 10}]

{{-10, 4.5653}, {-9, 4.72512}, {-8, 4.44479}, {-7, 4.18572}, {-6, 
  2.58606}, {-5, 3.00913}, {-4, 2.09279}, {-3, 3.71923}, {-2, 
  2.69165}, {-1, 1.55887}, {0, 4.69223}, {1, 3.07118}, {2, 
  5.43961}, {3, 3.50293}, {4, 2.9189}, {5, 3.70878}, {6, 3.57311}, {7,
   3.41627}, {8, 2.79404}, {9, 3.04386}, {10, 2.47882}}

You may ListPlot this if you like: my eyes can't easily see the underlying regularity in all of the noise. Nevertheless, Fit figures it out pretty well:

Fit[data, {1, UnitStep[x], x}, x]

2.10015 - 0.228675 x + 2.55593 UnitStep[x]

It works well here because my data represents a sum of known functions with unknown coefficients, and I fed those known functions to Fit as a basis. More phenomenological uses, like fitting a polynomial or Fourier series to data whose genesis is different, are trickier. High-order polynomials represented as power series (1,x^1,x^2,x^3,...) are especially tricky, as the results are generally numerically unstable.

"Function" in Mathematica has several different meanings. There are term rewriting formulae triggered by tags, like:

f[x_]:=x^2

Then there are "pure functions" made by wrapping Function around expressions built around Slot objects.

Function[Slot[1]^2]

For Fit, the "functions" are simply expressions containing zero or more undefined symbols that act as free variables. A "function" that lacks the free variable symbol(s) serves as a constant in the fit.