Write \begin{align} A &= \sum_{k=1}^r C_k E_k \in \mathbb{R}^{n \times n} \end{align} where $C_k$ are diagonal matrices and $E_k$ are circulant matrices for all $k \in \{1,2,\dots,r\}$. Given $x \in \mathbb{R}^n$, we can calculate $Ax$ and $A^Tx$ using the Fourier transform, that is \begin{align} Ax &= \sum_{k=1}^r C_k E_k x = \sum_{k=1}^r C_k \frac{1}{n} F^{-1} D_k F x\\ A^T x &= A^* x = \sum_{k=1}^r E_k^* C_k x = \sum_{k=1}^r \frac{1}{n} F^{-1} D_k^* F C_k x, \end{align} where $F$ is the DFT matrix, and $D_k$ is a diagonal matrix with diagonal elements given by the fourier transform of the first column of $E_k$, for all $k$.
My question is: can we calculate the product $A^T A$ using a similar strategy in the frequency domain? I tried to diagonalize $A^T A$ by isolating the product \begin{align} \frac{1}{n} F A^T A F^{-1} &= \sum_{k=1}^r \sum_{l=1}^r D_k^* F C_k C_l \frac{1}{n} F^{-1} D_l, \end{align} but I don't think it is a diagonal matrix, and so I don't think $A^T A$ is circulant. If it is not possible to calculate $A^TA$ in the frequency domain directly, can we at least approximate it using the operations $Ax$ and $A^Tx$ somehow?
Disclaimer: To make things a little bit more clear, in my application I don't have direct access to the matrix $A$ itself, only the first columns of each $E_k$, which are the convolution kernels. Also, I actually simplified the problem to 1d convolutions to write this question, but I'm working with 2d convolutions of vectorized images $x$, so each $E_k$ is doubly block circulant. The explicit construction of $A$ is very complicated and costly, that's why I'm avoiding it.
The matrix $A^T A$ is indeed circulant, since it is the product of two circulant matrices and circulant matrices form a ring. If $A$ has first row $(a_0,a_1,\dots, a_{n-1})$, then $A_T$ has first row $(a_0,a_{n-1},\dots,a_1)$. Moreover, the eigenvalues of $A$ are given by $\lambda_k=\sum_{i=0}^{n-1}a_i \omega^{ki}$ for $k=0,\dots, n-1$ and the corresponding eigenvalues for $A^T$ are given by $\tilde{\lambda_k}=\sum_{i=0}^{n-1}a_{n-i} \omega^{ki}=\lambda_{n-k}$.
So $A$ is diagonalized as $A=F D F^{-1}$, with $D=\text{diag}(\lambda_0,\dots, \lambda_{n-1})$, and $A^T$ is diagonalized as $A^T=F\tilde D F^{-1}$ where $\tilde D=\text{diag}(\lambda_{0},\lambda_{n-1},\dots, \lambda_{1})$. It follows that
$$A^T A=F \tilde D D F^{-1}$$ where the $\tilde D D=\text{diag}(\lambda_0^2,\lambda_1\lambda_{n-1},\dots, \lambda_1\lambda_{n-1})$. Note also that $\lambda_{n-i} =\overline{\lambda_i}$ so the eigenvalues of $A^TA$ are real, as expected from a symmetric matrix.
The above works only if $A$ is circulant. In general, if $E_k$ has first row $(e_0^k,\dots,e_{n-1}^k)$ then we can write $E_k=\sum_{i=0}^{n-1} e_i^kP^i$ in terms of the circulant matrix $P$ with first row $(0,1,0,\dots,0)$. Note that $E_k^T=\sum_{i=0}^{n-1} e_i^kP^{-i}$ and that $$P^j \text{diag}(c_0,c_1,\dots, c_{n-1})= \text{diag}(\sigma^j(c_0,c_1,\dots, c_{n-1})) P^j$$ where $\sigma$ shifts the elements in $(c_0,\dots, c_{n-1})$ one step to the left. Then we get \begin{align} A^TA &=\sum_{l,k=1}^r E_k^T C_k C_l E_l\\ &=\sum_{l,k=1}^r (\sum_{i=0}^{n-1} e_i^kP^{-i}) C_k C_l(\sum_{j=0}^{n-1} e_j^lP^j)\\ &=\sum_{i,j=0}^{n-1}\sum_{l,k=1}^r e_i^k e_j^l \sigma^{-i}(C_k C_l)P^{j-i}. \end{align} This shows that $A^TA$ is again of the form $\sum_i \tilde C_i \tilde E_i$ for diagonal matrices $\tilde C_i$ and circulant matrices $\tilde E_i$. But I don't see a way to use the same strategy as above if $A$ is not circulant.