Reputation: 1375
Let A
be a block circulant matrix with circulant blocks (i.e a BCCB matrix):
A = [1 2 3 4
2 1 4 3
3 4 1 2
4 3 2 1]
that is:
A = [C1 C2
C2 C1]
where each block (C1
, C2
) is a circulant matrix. I've read (see here) that BCCB can be diagonalized by following the equation: A =F*·D·F
where F
is the 2-D discrete Fourier transform matrix, F*
is the conjugate of F
, and D
is a diagonal matrix whose entries are the eigenvalues of A
.
In MATLAB I use this code:
(conj(dftmtx(4))/16*(fft2(A))*dftmtx(4))
but the result is:
[1 4 3 2
2 3 4 1
3 2 1 4
4 1 2 3]
Here the second and the fourth columns of A
are switched. Where is the error?
Upvotes: 5
Views: 9345
Reputation: 32930
Your source is a bit misleading. Diagonalizing a BCCB matrix with DFT is done as follows:
A = (FM⊗FN)*D(FM⊗FN)
where FN is the N-point DFT matrix, M is the number of Cj blocks and N is the size of each individual block (in your example M=2 and N=2). The "⊗" symbol denotes the tensor product.
Also note that F* = conj(F)T
(F* is called the complex conjugate transpose matrix). In MATLAB it translates to F'
instead of conj(F)
. Coincidentally, the DFT matrix F
is symmetric, which means that F* = conj(F)
is also true.
I'm not sure what you are trying to compute, but here's how the diagonalization of A
is done in MATLAB:
M = 2; N = 2;
FF = kron(dftmtx(M), dftmtx(N)); %// Tensor product
D = FF' * A * FF / size(A, 1); %// ' is the conjugate transpose operator
which yields:
D =
10 0 0 0
0 -2 0 0
0 0 -4 0
0 0 0 0
To diagonalize A
using only 2-D FFT operations, you can do this instead:
c = reshape(A(:, 1), N, []); %// First column of each block
X = fft2(c);
D = diag(X(:));
or in a one-liner:
D = diag(reshape(fft2(reshape(A(:, 1), N, [])), [], 1));
All of these produce the same diagonal matrix D
.
Hope this clarifies things for you!
Upvotes: 6