help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Complex matrix multiplication fails


From: Børge Strand-Bergesen
Subject: Complex matrix multiplication fails
Date: Wed, 12 May 2010 14:21:52 +0200

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


reply via email to

[Prev in Thread] Current Thread [Next in Thread]