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:
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.
�����������������������������������������������������������������������
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.
�����������������������������������������������������������������������
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.
�����������������������������������������������������������������������
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 :)