Circular Convolution Matrix of $ {H}^{H} {H} $

Question

We all know that Discrete Fourier Transform (DFT) corresponds to circular (not linear) convolution. That is to say, if $x(n),h(n)$ and $y(n)$ is the original signal, the filter and output signal in domain and $X(\omega),H(\omega),Y(\omega)$ are their counterpart in Fourier domain, then if we want $Y(\omega) = H(\omega)X(\omega)$ we must have $y(n) = h(n) \circledast x(n)$, where $\circledast$ is circular convolution (modulo the length of x).

If $h$ is the linear operator that stands for circular convolution with $h(n)$, then I think $h^T$, the transpose of $h$, is a circular convolution plus a circular shift. My question is, is $h^Th$ also a circular convolution of some kernel?

The motivation for thinking about this is to look for a fast method to compute $(I+h^Th)^{-1}$ through FFT. I know it seems common to introduce FFT to compute this, but how exactly when h is a convolution kernel?

Thanks in advance, Bob

Answer 1

If $ H $ is a matrix form of Circular Convolution then it is a Circulant Matrix.
Being a Circulant Matrix means it can be diagonalized by the Fourier Matrix $ {F} $:

$$ H = {F}^{H} D F $$

Where the matrix $ D $ id a Diagonal Matrix with the Fourier Coefficients of $ \mathcal{F} \left( h \right) $ on it main diagonal.
Also pay attention that we use the Conjugate Transpose $ {\left( \cdot \right)}^{H} $ operator as $ F $ and $ D $ are built by complex numbers (The space is over the complex numbers, hence the Adjoint Operator is the Complex Conjugate).

Now, it is easy to see that:

$$ {H}^{H} = {F}^{H} {D}^{H} F $$

Which is again an operator diagonalized with the Fourier Matrix hence it is a Circulant Matrix -> Circular Operation.
Basically it is a multiplication in Fourier Domain (DFT) by the conjugate of the DFT of the vector $ h $. By the DFT properties it means it is the flipped version of $ h $.

By the way:

$$ {H}^{H} H = {F}^{H} {D}^{H} F {F}^{H} D F = {F}^{H} {D}^{H} D F $$

Now, pay attention to $ {D}^{H} D $. It is a diagonal matrix multiplied by its conjugate. Namely it is, in the Fourier Domain, the Squared Magnitude.
In time domain it means it can be done by circular convolution with the correlation function of $ h $ (Or just do circular convolution in one direction and then in the other which is basically convolution with the flipped version).

Answer 2

If A is a "full" linear convolution operation (|signal|+|kernel|-1), A^T is a "valid" cross-correlation (shorter output, |signal| - |kernel| + 1). So I suppose

I suppose your h^T is a circular cross-correlation.

I am guessing you represented your convolution kernel h to a matrix H with Toeplitz structure and you are solving for g = (I+H^T H)^(-1) f thru iterative methods like conjugate-gradient which involves a lot of matrix-vector calculation. In particular you need to be able to calculate this H^T H v fast for any vector v ? Then you can exploit circular Toeplitz structure for fast matrix-vector calculation. Wikipedia has some reference on that.