CASTalk.com

vivi2003 · Guest

This is a program I found in the book named: DSP using Matlab v4.0. I thin
it's quite popular since I can easily find in on the Internet using google
Many website provide its source code

It's about Circular convolution:
----------------------------------------------------
function y = circonvt(x1,x2,N)
if length(x1) > N
error('N must be >= the length of x1')
end
if length(x2) > N
error('N must be >= the length of x2')
end
x1=[x1 zeros(1,N-length(x1))];
x2=[x2 zeros(1,N-length(x2))];

m = [0:1:N-1];
x2 = x2(mod(-m,N)+1);
H = zeros(N,N);
for n = 1:1:N
H(n,:) = cirshftt(x2,n-1,N);
end
y = x1*H';
-----------------------------------------------------

I think that this program has a bug since it give different results with:
x1=rand(1,N)+j*rand(1,N);
x2=rand(1,N)+j*rand(1,N);

-> y1=circonvt(x1,x2,N) <> y2=circonvt(x2,x1,N) ???

meanwhile circular convolution is commutative so we expect that the resul
should be identical.

Any one help me plz :)

Jerry Avins · Guest

vivi2003 wrote:

Jerry Avins · Guest

Praveen wrote:

Praveen · Guest

Hey rand is a random function, so obviously you are bound to get
different results!!
Relax. You are doing just fine!!!

Praveen

David Tweed · Guest

vivi2003 wrote:

Jerry Avins · Guest

David Tweed wrote:

vivi2003 · Guest

Thank you very much for the attention. I wrote directly to the author, Pro
V.Ingle and here is his question. Should any one face the problem as I do
plz read this :)

--------------------------
Dear Prof. NGUYEN:

Thanks for the feedback. Yes, there appears to be a bug. The function
was
originally developed in mid 1980's when matrix transpose H' was a
transpose
operator. Sometime later along the way, Matlab changed the matrix
transpose
operator and made it the conjugate transpose operator. This has caused
problems
for some of my old functions which I thought I has corrected. The
circonvt
function calls for the cirshift function which used H' for matrix
vector
multiplication. For real signals, the function worked but I missed the
circonvt
function for complex signals. A simple fix is to change the last line

y = x1*H';

to

y = x1*conj(H');

The complete function is
----------------------------------------------------------------------------------
function y = circonvt(x1,x2,N)
% N-point circular convolution between x1 and x2: (time-domain)
% -------------------------------------------------------------
% [y] = circonvt(x1,x2,N)
% y = output sequence containing the circular convolution
% x1 = input sequence of length N1 <= N
% x2 = input sequence of length N2 <= N
% N = size of circular buffer
% Method: y(n) = sum (x1(m)*x2((n-m) mod N))

% Check for length of x1
if length(x1) > N
error('N must be >= the length of x1')
end

% Check for length of x2
if length(x2) > N
error('N must be >= the length of x2')
end

x1=[x1 zeros(1,N-length(x1))];
x2=[x2 zeros(1,N-length(x2))];

m = [0:1:N-1];
x2 = x2(mod(-m,N)+1);
H = zeros(N,N);

for n = 1:1:N
H(n,:) = cirshftt(x2,n-1,N);
end

y = x1*conj(H');
---------------------------------------------------------------------------------

Verification:

	CASTalk.com Forum Index -> DSP	All times are GMT
Page 1 of 1