This approach uses "Number Walls".

Given a sequence of elements
$\,a:\mathbb{Z}\to F\,$ in a field $\,F.\,$ Define the infinite matrix
of $\,n\times n\,$ Hankel determinants

$$ W_{n,m} :=\det({a_{m-n+i+j}})_{i,j=0}^{n-1}. \tag{1} $$

According to a theorem on determinants by Jacobi, $\,W_{n,m}\,$ satisfies
the Number Wall identity (except at $n=m=0$)

$$ W_{n+1,m}W_{n-1,m} = W_{n,m+1}W_{n,m-1} - W_{n,m}W_{n,m} \tag{2} $$

if $\,a_n=0\,$ for all $\,n<0.\,$
The definition in equation $(1)$ implies

$$ W_{0,m} = 1,\quad W_{1,m} = a_{m-1}. \tag{3} $$

Define the sequence of functions

$$ p_n(k) := W_{k,n+k} = \det({a_{n+i+j}})_{i,j=0}^{k-1}. \tag{4} $$

Use this definition and equation $(3)$ to get

$$ p_n(0) = 1, \quad p_n(1) = a_n. \tag{5} $$

For simplicity, suppose $\,p_n\,$ is a polynomial of degree $\,n\,$ for
$\,n=0,1.\,$ Then equation $(5)$ implies

$$ p_0(x) = a_0 = 1, \quad p_1(x) = x(a_1-1)+1. \tag{6} $$

The definition in equation $(4)$ implies $\, p_0(k) = 1 = W_{k,k}\,$
which is a linear equation in $\,a_{2k}\,$ which can be solved for when
given values of all of the $\,a_0,a_1,\dots,a_{2k-1}.\,$

Similarly, $\,p_1(k) = k(a_1-1)+1 = W_{k,k+1}\,$ which is a linear
equation in $\,a_{2k+1}\,$ which can be solved for when given values of
all the $\,a_0,a_1,\dots,a_{2k}.\,$

Using this procedure, solve the equations successively to get (with $\,z:=a_1$)

$$ a_1 \!=\! z, \;\; a_2 \!=\! 1+z^2\!, \;\; a_3 \!=\!
2+2z+z^3\!, \;\; a_4 \!=\! 6+4z+3z^2+z^4,\;\;\dots. \tag{7} $$

An alternative method is to rewrite the Number Wall identity equation $(2)$ as

$$ p_{n+1}(k)\,p_{n-1}(k) = p_{n+1}(k-1)\,p_{n-1}(k+1) + p_n(k)\,p_n(k). \tag{8} $$

This implies the $\,p_n(k)\,$ recursion

$$ p_n(k) = (p_n(k-1)\,p_{n-2}(k+1) + p_ {n-1}^2(k))/p_{n-2}(k). \tag{9} $$

Now use equation $(5)$ to get $\,a_n = p_n(1).$

The nice polynomial expression for $\,a_n\,$ is

$$ a_n = \sum_{k=0}^n T_{n,k}\, z^k,\quad
\frac2{1+2x+\sqrt{1-4x}-2xy} = \sum_{n,k} T_{n,k}\, x^n y^k \tag{10} $$

where $\,T\,$ is the triangular array in
OEIS sequence A065600

Triangle T(n,k) giving number of Dyck paths of length 2n with exactly k hills (0 <= k <= n).

For $\,z=0\,$ we get OEIS sequence A000957 with
a different offset.

For $\,z=1\,$ we get $\,a_n=C_n.\,$

For $\,z=2\,$ we get $\,a_n=C_{n+1}.\,$

For $\,z=3\,$ we get $\,a_n=\binom{2n+1}{n}.\,$

For $\,z=4\,$ we get OEIS sequence A049027.

All of these sequences with $\,z>0\,$ are the rows of the infinite array in
OEIS sequence A076037 (with the first column of
all 1s removed).

Note that
$$ p_2(x) = \frac16(1+x)((6-5x+2x^2)-(4x-4x^2)z+(x+2x^2)z^2). \tag{11} $$
This seems to factor into linear factors only if $\,z=1,2,3.$
$$ p_3(x) = \frac1{180}(1+x)(2+x)(3+2x)F_3(x,z) \tag{12} $$ where
$\,F_3(x,z)\,$ is a cubic in $\,x\,$ and $\,z\,$ which also seems to
factor into linear factors only if $\,z=1,2,3.$
$$ p_4(x) = \frac1{75600}(1+x)(2+x)^2(3+x)(3+2x)(5+2x)F_4(x,z) \tag{13} $$ where
$\,F_4(x,z)\,$ is degree $4$ in $\,x\,$ and $\,z.\,$
It seems to factor into linear factors only if $\,z=1,2,3.$

This is consistent with the definitions of $\,f_n(x)\,$ and $\,g_n(x)\,$
given in the question.