There is a nice paper Enumeration of Finite Automata by Frank Harary and Ed Palmer which presents a formula $a(n,k,m)$ for the number of finite automata with $n$ states, $k$ input symbols and $m$ output symbols. It is stated in Corollary $3$ as \begin{align*} a(n,k,m)=\frac{1}{n!k!m!}\sum_{H_1} I(\alpha,\beta,\alpha)I(\alpha,\beta,\gamma) \end{align*} where the sum is over all permutations in $H_1=S_n^{S_n\times S_k}\times S_m^{S_n\times S_k}$ of the form $\{[(\alpha,\beta);\alpha^{-1}],[(\alpha,\beta);\gamma]\}$ with \begin{align*} I(\alpha,\beta,\gamma)=\prod_{p=1}^n\prod_{q=1}^k\left[\sum_{s|[p,q]}sj_s(\gamma)\right]^{j_p(\alpha)j_q(\beta)\langle p,q\rangle} \end{align*}
Here the authors denote with $[p,q]:=\operatorname{lcm}(p,q), \langle p,q\rangle:=\operatorname{gcd}(p,q)$ and $j_p(\alpha)$ denotes the number of cycles of the permutation $\alpha$ with length $p$.
The special case $n=k,m=1$ is already analysed and calculated for small values of $n$ in this MSE post, especially with focus to $n=4$.
Question: Maybe someone could provide some computations for the general case for small values of $n,k,m$?
In answering this question we refer to the algorithm at the MSE link which works for the generalized problem as well. The only difference is that the values that go into the slots of the array/table are pairs of states and output symbols, meaning when we transition from a certain column on an input symbol corresponding to a row we transition to the state (first element of the pair) and output the symbol (second element of the pair). The action on the slots is the simultaneous action of $\pi$ and $\tau$ on the rows and columns and we now have a permutation $\sigma$ which acts on the set of output symbols and the action on the values is the combined action of $\tau$ and $\sigma$ on the state / symbol pairs.
We get the following table for one output symbol.
For two output symbols we have
Three output symbols yield
Four symbols yield
Five yield
Finally we get for six symbols
The Maple code for this was as follows.
with(combinat); pet_cycleind_symm := proc(n) local p, s; option remember; if n=0 then return 1; fi; expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n)); end; pet_flatten_term := proc(varp) local terml, d, cf, v; terml := []; cf := varp; for v in indets(varp) do d := degree(varp, v); terml := [op(terml), seq(v, k=1..d)]; cf := cf/v^d; od; [cf, terml]; end; cycles_prod := proc(cyca, cycb) local ca, cb, lena, lenb, res, vlcm; res := 1; for ca in cyca do lena := op(1, ca); for cb in cycb do lenb := op(1, cb); vlcm := lcm(lena, lenb); res := res*a[vlcm]^(lena*lenb/vlcm); od; od; res; end; automaton := proc(N, M, K) option remember; local idx_slots, idx_cols, idx_syms, res, a, b, c, sim, flat_sim, sym, flat_sym, flat_a, flat_b, flat_c, cyc_a, cyc_b, len_a, len_b, p, q; if N > 1 then idx_slots := pet_cycleind_symm(N); else idx_slots := [a[1]]; fi; if M > 1 then idx_cols := pet_cycleind_symm(M); else idx_cols := [a[1]]; fi; if K > 1 then idx_syms := pet_cycleind_symm(K); else idx_syms := [a[1]]; fi; res := 0; for a in idx_slots do flat_a := pet_flatten_term(a); for b in idx_cols do flat_b := pet_flatten_term(b); sim := cycles_prod(flat_a[2], flat_b[2]); flat_sim := pet_flatten_term(sim); for c in idx_syms do flat_c := pet_flatten_term(c); sym := cycles_prod(flat_b[2], flat_c[2]); flat_sym := pet_flatten_term(sym); p := 1; for cyc_a in flat_sim[2] do len_a := op(1, cyc_a); q := 0; for cyc_b in flat_sym[2] do len_b := op(1, cyc_b); if len_a mod len_b = 0 then q := q + len_b; fi; od; p := p*q; od; res := res + p*flat_a[1]*flat_b[1]*flat_c[1]; od; od; od; res; end; output := proc(MXN, MXM, K) local data, N, M, fd, fname, width; data := table(); for N to MXN do data[N] := table(); for M to MXM do data[N][M] := automaton(M, N, K); od; od; fname := sprintf("automata-%d-%d-%d.txt", MXN, MXM, K); fd := fopen(fname, WRITE); for N to MXN do fprintf(fd, "|"); for M to MXM do width := nops(convert(data[MXN][M], base, 10)); fprintf(fd, "% *d|", width+1, data[N][M]); od; fprintf(fd, "\n"); od; fclose(fd); end;