Charcoal, 41 35 bytes

Ｗθ«Ｉ…θ¹≔ΦＥθ⁻κקθ⁰§↨Ｘ⊕Ｘ⁴Ｌθ⊖ＬθＸ²Ｌθλλθ

Try it online! Link is to verbose version of code. Takes the first n+1 coefficients as input. Explanation: Now a port of @Albert.Lang's current Python answer.

Ｗθ«

Repeat until all of the coefficients have been processed.

Ｉ…θ¹

Output the next coefficient.

≔ΦＥθ⁻κקθ⁰§↨Ｘ⊕Ｘ⁴Ｌθ⊖ＬθＸ²Ｌθλλθ

Divide by (x²+1)ⁿ (even without vectorised subtraction, most of the work here is performing the binomial calculation) then discard the first coefficient.

(Although vectorised subtraction can be emulated on the experimental zip branch, the byte count is the same, although there is a performance boost.)

Python, 89 bytes

-2 thanks to @Ted.

f=lambda p,n:p[:1]+(n*[c:=(4**n+1)**n]and f([d-(c:=c>>n)%2**n*p[0]for d in p[1:-1]],n-1))

Attempt This Online!

Previous Python, 103 bytes

f=lambda p,c=0:p and p[:1]+f([d-c%2*p[0]*math.comb(len(p)//2,(c:=c+1)//2)for d in p[1:-1]])
import math

Attempt This Online!

How?

Essentially, divides the polynomial \$p\$ by \$(x²+1)^n\$ keeps the quotient, replaces \$p\$ with the remainder, decrements \$n\$ and starts over.

Python3, 500 bytes

S=sorted
def M(*p):
 p,*q=p
 d=dict(p)
 for Q in q:
  D={}
  for i,a in d.items():
   for I,A in Q:D[i+I]=D.get(i+I,0)+a*A
  d=D
 return d.items()
def f(p):
 Q,s=[(V:=[[len(p)//2,0]],V[0][0])],[V]
 for q,n in Q:
  if S(M([(len(p)//2,1)],[j for a,b in q for j in M([(0,b)],*[[(1,1),(-1,1)]for _ in range(a)])]))==S(p):return q
  U=[(q+[[n-1,0]],n-1),(q,n-1)]*(n>0)
  for i in range(len(q)):
   for v in[1,-1]:u=eval(str(q));u[i][1]+=v;U+=[(u,n)]
  for u,n in U:
   if S(u)not in s:Q+=[(u,n)];s+=[S(u)]

Try it online!

Wolfram Language (Mathematica), 38 bytes

Simplify[#/x^#2/.x->(x+√(x^2-4))/2]&

Try all test cases online! An unnamed function taking a polynomial of degree \$2n\$ as its first argument, formatted like x^4+2x^3+3x^2+2x+1, and the semidegree \$n\$ as its second argument. (Note that √ is a 3-byte character.) The output is given by a polynomial, sometimes in factored form like (1+x)^2 instead of 1+2x+x^2; if that's not sufficient, prepend Expand@ to the code for an additional 7 bytes.

#/x^#2 gives the polynomial divided by the appropriate power of \$x\$, which is the thing (call it \$h(x)\$) that we want to equal \$q(x+\frac1x)\$. But \$(x+\sqrt{x^2-4})/2\$ is the inverse function to \$x+\frac1x\$, so \$h(x)=q(x+\frac1x)\$ is equivalent to \$q(x)=h\bigl((x+\sqrt{x^2-4})/2\bigr)\$, which is what the code computes. Simplify is needed to multiply out all the powers of \$\sqrt{x^2-4}\$ and get them to collapse into regular polynomials.

(For the mathies: \$x+\frac1x\$ has a second inverse function \$(x-\sqrt{x^2-4})/2\$ on a different domain, but it yields the same answer, as we could see either by analytic continuation mumble mumble or Galois theory mumble mumble.)

x86-64 machine code, 26 bytes

AD 31 C9 EB 08 FF C1 F7 F9 29 44 CE FC F7 EF FF CF 39 F9 7C F0 01 CF 79 E7 C3

Try it online!

Following the standard calling convention for Unix-like systems (from the System V AMD64 ABI), this takes \$n\$ in EDI and the address of an array of signed 32-bit integers in RSI (the first \$n+1\$ coefficients), and modifies the array in place.

This works by subtracting \$\binom{n}{i}a_0\$ from \$a_{2i}\$ for \$0 < i \leq \frac{n}{2}\$ (those are the coefficients from expanding \$a_0(x+\frac{1}{x})^n\$), and repeating the process on the subarray without the first element. Each value of \$\binom{n}{i}a_0\$ is calculated from the previous one, by multiplying by \$(n - i + 1)\$ and dividing by \$i\$.

In assembly:

f:  lodsd           # Load the first value into EAX, advancing RSI.
    xor ecx, ecx    # Set ECX to 0.
    jmp a           # Jump into the middle of the loop.
r:  inc ecx                     # Add 1 to ECX.
    idiv ecx                    # Divide the product by ECX (quotient in EAX).
    sub [rsi + 8*rcx - 4], eax  # Subtract EAX from the appropriate value in memory.
a:  imul edi        # Multiply EAX by EDI (initially the degree).
    dec edi         # Subtract 1 from EDI.
    cmp ecx, edi    # Compare ECX and EDI.
    jl r            # Jump if ECX is lesser.
    add edi, ecx    # Add ECX to EDI. EDI is 1 less than it was initially.
    jns f           # Jump back to the start if that value is nonnegative.
    ret             # Return.

Python, 128 bytes

lambda p,n:map(round,numpy.polyfit(r:=range(3,2*n+4),[sum(c*((a+(a**2-4)**.5)/2)**(i-n)for i,c in p)for a in r],n))
import numpy

Attempt This Online!

Explanation

We want to find \$q(x)\$. From \$p(x) = x^nq(x+1/x)\$, we know that \$q(x + 1/x) = p(x)/x^n\$

We can calculate \$p(x)/x^n\$ at a given \$x\$ by shifting the coefficients of \$p(x)\$, and then substituting in a value for \$x\$, as below.

sum([c*(<x>)**(i-n)for i,c in p]

However we actually want to transform \$x\$ in order to invert the \$x+1/x\$ which has been applied to \$q\$. One suitable function which does this is \$\frac{x+\sqrt{x^2-4}}{2}\$. Note that the domain for this function excludes \$[-2,2]\$ as \$x^2-4\$ must be positive.

So we apply the transformation (a+(a**2-4)**.5)/2 (a is representing a value of x), and sample values from 3 to (order of p) + 4 (excluding the final value due to how range works).

We can then use these values to fit a polynomial of the correct order to the data points we've created. This is \$q(x)\$.