Reply by vivi2003 December 4, 20052005-12-04
Thank you very much for the attention. I wrote directly to the author, Prof
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:

>> N=10; >> x1=rand(1,N)+j*rand(1,N); >> x2=rand(1,N)+j*rand(1,N); >> circonvt(x1,x2,N)
ans = Columns 1 through 4 -1.7958 + 4.3439i -1.5197 + 5.3775i -1.9221 + 4.4140i -2.3498 + 5.1267i Columns 5 through 8 -2.2629 + 5.0657i -2.0736 + 4.8734i -1.6936 + 4.6140i -1.3524 + 5.1862i Columns 9 through 10 -2.5356 + 4.9775i -2.2091 + 5.2821i
>> circonvt(x2,x1,N)
ans = Columns 1 through 4 -1.7958 + 4.3439i -1.5197 + 5.3775i -1.9221 + 4.4140i -2.3498 + 5.1267i Columns 5 through 8 -2.2629 + 5.0657i -2.0736 + 4.8734i -1.6936 + 4.6140i -1.3524 + 5.1862i Columns 9 through 10 -2.5356 + 4.9775i -2.2091 + 5.2821i
>> ifft(fft(x1,N).*fft(x2,N))
ans = Columns 1 through 4 -1.7958 + 4.3439i -1.5197 + 5.3775i -1.9221 + 4.4140i -2.3498 + 5.1267i Columns 5 through 8 -2.2629 + 5.0657i -2.0736 + 4.8734i -1.6936 + 4.6140i -1.3524 + 5.1862i Columns 9 through 10 -2.5356 + 4.9775i -2.2091 + 5.2821i
>> max(abs(circonvt(x1,x2,N)-circonvt(x2,x1,N)))
ans = 2.8087e-015 Thanks again for pointing this out. Regards, Prof. Vinay K. Ingle
Reply by Jerry Avins December 3, 20052005-12-03
David Tweed wrote:
> vivi2003 wrote: > >>This is a program I found in the book named: DSP using Matlab v4.0. >>I think it's quite popular since I can easily find in on the Internet >>using google. Many website provide its source code > > > Ignore Jerry and Praveen; clearly they have not read your message > well enough to understand what you're doing. > > I looked at one version of the source code for circonvt() I found > on the web (it matches yours, but has some additional comments), > but I have to confess that I don't know Matlab well enough to > explain what the problem is. > > I strongly suspect it has something to do with the odd way in which > x2 gets sort-of reversed in-place: > > m = [0:1:N-1]; > x2 = x2(mod(-m,N)+1); > > Note that the first element of x2 doesn't move. Could the second > line be written as "x2 = x2(N-m);" instead?
David, You're right. I didn't understand because I didn't read carefully enough. My question to Praveen stands. Jerry -- Engineering is the art of making what you want from things you can get. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
Reply by David Tweed December 3, 20052005-12-03
vivi2003 wrote:
> This is a program I found in the book named: DSP using Matlab v4.0. > I think it's quite popular since I can easily find in on the Internet > using google. Many website provide its source code
Ignore Jerry and Praveen; clearly they have not read your message well enough to understand what you're doing. I looked at one version of the source code for circonvt() I found on the web (it matches yours, but has some additional comments), but I have to confess that I don't know Matlab well enough to explain what the problem is. I strongly suspect it has something to do with the odd way in which x2 gets sort-of reversed in-place: m = [0:1:N-1]; x2 = x2(mod(-m,N)+1); Note that the first element of x2 doesn't move. Could the second line be written as "x2 = x2(N-m);" instead? -- Dave Tweed
Reply by Jerry Avins December 3, 20052005-12-03
Praveen wrote:
> Hey rand is a random function, so obviously you are bound to get > different results!! > Relax. You are doing just fine!!! > > Praveen
Praveen, I must have missed something. What did you intend your message to add to mine? Jerry -- Engineering is the art of making what you want from things you can get. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
Reply by Praveen December 3, 20052005-12-03
Hey rand is a random function, so obviously you are bound to get
different results!!
Relax.  You are doing just fine!!!

Praveen

Reply by Jerry Avins December 3, 20052005-12-03
vivi2003 wrote:
> This is a program I found in the book named: DSP using Matlab v4.0. I think > 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 result > should be identical. > > Any one help me plz :)
Plz yourself! Don't you get a new result on every call to rand()? Why do you expect the same result with different arguments? Jerry -- Engineering is the art of making what you want from things you can get. &#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
Reply by vivi2003 December 3, 20052005-12-03
This is a program I found in the book named: DSP using Matlab v4.0. I think
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 result
should be identical. 

Any one help me plz :)