Re: Complex matrix multiplication fails

[Børge Strand-Bergesen]

Just a side note: The matrix is symmetrical, so that multiplication by
CEM(:,k+1)) is the same as multiplication by CEM(k+1,:)). The code
behaves the same way regardless of swapping the indexes.

Borge


On Wed, May 12, 2010 at 14:21, Børge Strand-Bergesen
<address@hidden> wrote:
> Hi,
>
> I'm playing with some Discrete Fourier Transform code that I'd like to
> implement as a complex
> matrix multiplication. When I write out the complete sequence or apply
> the matrix row-by-row,
> everything looks fine. But when I perform the complex matrix
> multiplication, the results do not
> match. I have performed a check with a real matrix for the same input,
> where the results match.
>
> Do you see any errors in my code? Are there any known bugs here?
> Curiously, if I replace
> CEM=WN.^CEM; by CEM=WM.^-CEM; the complex matrix multiplication
> matches fft() but the matrix
> applied manually does not.
>
> I'm using Octave 3.2.2 natively on Win7/32.
>
>
> Thanks,
> Borge
>
>
> % Matrix-based DFT
> % n and k are counted from 0 to N-1. An offset of 1 is added for each
> matrix index.
>
> % Our test filter in the time domain
> xn=[1 1 0.5 0 0 0 0 0];
> N=length(xn);
> WN=exp(-j*2*pi/N);      % Complex exp. shorthand
>
> % Written-out DFT === OK ===
> for k=0:N-1, Xk(k+1)=0; for n=0:N-1,
> Xk(k+1)=Xk(k+1)+xn(n+1)*WN.^(n*k); end; end;
> max(abs([imag(Xk-fft(xn)) real(Xk-fft(xn))]))
> % ans = 1.3323e-015
>
> % Alternatively === OK ===
> for k=0:N-1, Xk(k+1)=sum(xn.*WN.^([0:N-1]*k)); end;
> max(abs([imag(Xk-fft(xn)) real(Xk-fft(xn))]))
> % ans = 1.3323e-015
>
> % Matrix-based DFT, complex exponential matrix
> CEM=[0:N-1]'; for n=1:N-1, CEM(:,n+1)=n*CEM(:,1); end; CEM(:,1)=0*[0:N-1];
> % === CRUCUAL DIFFERENCE ===
> CEM=WN.^CEM; % Matches WNN=exp(-j*2*pi/N*WNN);
>
> % Matrix applied manually === OK with CEM=WN.^CEM; FAILS with CEM=WN.^-CEM; 
> ===
> for k=0:N-1, Xk(k+1)=sum(xn.*CEM(:,k+1)); end;
> max(abs([imag(Xk-fft(xn)) real(Xk-fft(xn))]))
> % ans = 1.3323e-015
>
> % Matrix multiplication === FAILS with CEM=WN.^CEM; OK with CEM=WN.^-CEM; ===
> Xk = CEM * xn'; Xk=Xk';
> max(abs([imag(Xk-fft(xn)) real(Xk-fft(xn))]))
> % ans =  2.4142
>
>
> % Repeated with real-value matrix
>
> CEM=[0:N-1]'; for n=1:N-1, CEM(:,n+1)=n*CEM(:,1); end; CEM(:,1)=0*[0:N-1];
>
> % Matrix applied manually
> for k=0:N-1, Xk1(k+1)=sum(xn.*CEM(:,k+1)); end;
>
> % Matrix multiplication
> Xk2 = CEM * xn'; Xk2=Xk2';
> max(abs(Xk1-Xk2))
> % ans = 0
>