I've figured how to do this with MAGMA. I'm going to document it, in case anyone has the same problem in the future:
1) Define a parameter field:
R$<$a1,...,an$>$:=FieldOfFractions(PolynomialRing(RationalField(),n));
where a1,...,an are the parameters and n is their number.
2) Define the free algebra:
F$<$x1,...xn,y1,...,yn$>$:=FreeAlgebra(R,2n);
where x1,...,xn are the variables and y1,...,yn are going to be their images through the required morphism (i.e., either y1=(x1)' or y1=(x1)^*, and so on).
3a) Involutive automorphism:
This case is straightforward, because MAGMA implements homomorphisms defined on generators and extended by linearity. Define the involutive automorphism via
h := hom$<$F -> F | y1,...,yn,x1,...,xn$>$;
3b) Involution:
Since antihomomorphisms are not implemented in MAGMA, we have to code the full involution. To do so:
3b1) Write a function taking one parameter (a nonconmmutative polynomial):
involution:=function(poly):
...
end function;
3b2) Do a loop, iterating every term of the polynomial:
m:=Terms(poly);
for term in m do
...
end for;
3b3) Extract the coefficient of the term and read the variables of the term backwards with another loop:
coef:=Coefficients(n)[1];
for i:=Length(term) to 1 by -1 do
...
end for;
3b4) Control the image of every variable with a case instruction:
case term[i]:
when x1: res:=y1;
...
end case;
3b5) After the case, keep multiplying the partial results in order to get the involuted monomial (after all iterations of the second loop).
3b6) At the end of the first loop, multiply the involuted monomial by the stored coefficient and keep adding the involuted terms in order to get the involuted polynomial.