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 :)
help me find the bug in this program please :-??
Started by ●December 3, 2005
Reply by ●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 ●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 ●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!!! > > PraveenPraveen, 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 ●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 codeIgnore 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 ●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 ●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