Reply by Rune Allnor April 18, 20092009-04-18
On 17 Apr, 23:18, itakatz <itak...@gmail.com> wrote:
> On Apr 17, 1:42&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote: > > > > > > > On 17 Apr, 09:54, itakatz <itak...@gmail.com> wrote: > > > > > 4) &#4294967295;Compute the DFT of the zero-padded frame > > > > 5) &#4294967295;Compute the periodogram as the squared magnitude > > > > &#4294967295; &#4294967295; of the spectrum coefficients > > > > 6) &#4294967295;Compute the IDFT of the periodogram > > > > > Rune > > > > A little bit off-topic: > > > > Since multiplication in frequency domain is the same as (circular) > > > convolving in the time domain, > > > Convolution by DFT is circular, which is why you need > > to zero-pad before computing the DFT. The convolution > > of two N-length sequences gives as result a sequence > > with at most 2N-1 non-zero elements. Zero-padding to > > at least 2N-1 total length allocates the space needed > > to avoid wrap-around effects. Depending on what > > implementation of the FFT you have available, you might > > zero-pad further, to the next power of 2. > > > > shouldn't one of the segments be time- > > > reversed before zero-padding and DFTing? > > > Time-reversing in time domain is not sufficient to > > handle wrap-around effects from circular convolution, > > since you don't change the length of the input sequence. > > > With pen and paper your trick might work, but in the > > computer you need to allocate memory and find a map > > from the 'true' index space k = -N,...,N to the RAM > > address space, BaseAddress + {0,1,...2N}. > > > Essentially, you need to do some book-keeping to find > > out if the rxx[0] coefficient turns out in the index 0 > > position of the buffer, or in the index N position. > > > > in that way the auto- > > > correlation is re-written as a convolution, and then it is ok to > > > compute it in the freq domain. Am I right? > > > In time domain, reversing the direction of a sequence x[n] > > is written as x[-n]. It's a standard excercise in most DSP > > books to show that if > > > x[n] <-> X[k] > > > is a DFT pair, then > > > x[-n] <-> X'[k] > > > where X' is the complex conjugate of X, is also a DFT pair. > > Well, this is correct at least when x[n] is real-valued, > > it might not be for complex-valued x[n]. > > > Anyway, convolving x[n] by x[-n] in time domain, which > > is a common estimator for the autocorrelation function > > (save some scaling coefficients), can alternatively be > > computed in frequency domain as a product of the DFTs: > > > Sxx[k] = X[k]X'[k] = |X[k]|^2. > > > Hence my algorithm works, once the relevant scaling > > coefficients have been accounted for. > > > Rune > > Thanks. > > In the meantime I did the simple math and I indeed got that the DFT of > the time reversed signal is the c.c. of the DFT of the original > signal, up to a linear phase, which accounts for the fact that my time > reversed signal was x[N-n-1] instead of x[-n].
The practical implication would probably be that you find the rxx[0] element of the autocorrelation in array element with index N, instead of the element with index 0. Not a big deal, but a detail to keep track of. Rune
Reply by itakatz April 17, 20092009-04-17
On Apr 17, 1:42&#4294967295;pm, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 17 Apr, 09:54, itakatz <itak...@gmail.com> wrote: > > > > 4) &#4294967295;Compute the DFT of the zero-padded frame > > > 5) &#4294967295;Compute the periodogram as the squared magnitude > > > &#4294967295; &#4294967295; of the spectrum coefficients > > > 6) &#4294967295;Compute the IDFT of the periodogram > > > > Rune > > > A little bit off-topic: > > > Since multiplication in frequency domain is the same as (circular) > > convolving in the time domain, > > Convolution by DFT is circular, which is why you need > to zero-pad before computing the DFT. The convolution > of two N-length sequences gives as result a sequence > with at most 2N-1 non-zero elements. Zero-padding to > at least 2N-1 total length allocates the space needed > to avoid wrap-around effects. Depending on what > implementation of the FFT you have available, you might > zero-pad further, to the next power of 2. > > > shouldn't one of the segments be time- > > reversed before zero-padding and DFTing? > > Time-reversing in time domain is not sufficient to > handle wrap-around effects from circular convolution, > since you don't change the length of the input sequence. > > With pen and paper your trick might work, but in the > computer you need to allocate memory and find a map > from the 'true' index space k = -N,...,N to the RAM > address space, BaseAddress + {0,1,...2N}. > > Essentially, you need to do some book-keeping to find > out if the rxx[0] coefficient turns out in the index 0 > position of the buffer, or in the index N position. > > > in that way the auto- > > correlation is re-written as a convolution, and then it is ok to > > compute it in the freq domain. Am I right? > > In time domain, reversing the direction of a sequence x[n] > is written as x[-n]. It's a standard excercise in most DSP > books to show that if > > x[n] <-> X[k] > > is a DFT pair, then > > x[-n] <-> X'[k] > > where X' is the complex conjugate of X, is also a DFT pair. > Well, this is correct at least when x[n] is real-valued, > it might not be for complex-valued x[n]. > > Anyway, convolving x[n] by x[-n] in time domain, which > is a common estimator for the autocorrelation function > (save some scaling coefficients), can alternatively be > computed in frequency domain as a product of the DFTs: > > Sxx[k] = X[k]X'[k] = |X[k]|^2. > > Hence my algorithm works, once the relevant scaling > coefficients have been accounted for. > > Rune
Thanks. In the meantime I did the simple math and I indeed got that the DFT of the time reversed signal is the c.c. of the DFT of the original signal, up to a linear phase, which accounts for the fact that my time reversed signal was x[N-n-1] instead of x[-n].
Reply by Rune Allnor April 17, 20092009-04-17
On 17 Apr, 16:11, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Apr 16, 1:39&#4294967295;am, "sasuke" <a.ssjg...@gmail.com> wrote: > > > I am using MATLAB to generate close to 10^8 samples of a Rayleigh fading > > channel using Jakes method. Once I generate these samples, I have to then > > create a 20 x 20 autocorrelation matrix. I am using MATLAB to do this. > > MATLAB works perfectly for 10^6 samples, but once the number of samples > > increases slightly above 50 * 10^6, MATLAB throws an "OUT of Memory" error. > > I tried some techniques given inwww.mathworks.comtoovercome this > > problem, but they don't seem to help once the number of samples goes above > > 10^7. > > You must be using a machine from several years ago, or are using an > algorithm that requires lots of extra memory. &#4294967295;50*10^6 double- > precision samples is under 400MB of memory.
Once upon a time matlab by default used algorithms based on verbatim implementations of matrix algebra expressions. In recent times that's not quite as much of a default, for canned functions supplied with matlab or its toolboxes, but this 'vectorization' nonsense certainly tends to lead naive users in that direction. With that in mind, I would expect 'experienced' naive users [*] to implement the computation of an K x K autocorrelation matrix from an N-length vector as something like 1) Build a K x N'ish data matrix X from the input data (maybe adding or subtracting end effects, depending on what estimator one implements) 2) Do the final computation as Rxx = X*X'; Which means one needs at least one, maybe two, copies of an K x N'ish matrix in memory to do the computations. With the numbers the OP mentions, that's 20 * 1e8 = 2e9 elements, or 1.6e10 bytes per copy of the matrix. With such numbers there's no wonder he runs out of memory. I don't have access to the signal processing toolbox, so I can't check it out myself. Does anybody know what algorithms matlab SPT functions actually use to compute the autocorrelation matrices? Rune [*] Users who know enough matlab to know about 'vectorization', but not enough about computers or programming to understand what they actually are doing.
Reply by Steven G. Johnson April 17, 20092009-04-17
On Apr 16, 1:39&#4294967295;am, "sasuke" <a.ssjg...@gmail.com> wrote:
> I am using MATLAB to generate close to 10^8 samples of a Rayleigh fading > channel using Jakes method. Once I generate these samples, I have to then > create a 20 x 20 autocorrelation matrix. I am using MATLAB to do this. > MATLAB works perfectly for 10^6 samples, but once the number of samples > increases slightly above 50 * 10^6, MATLAB throws an "OUT of Memory" error. > I tried some techniques given inwww.mathworks.comto overcome this > problem, but they don't seem to help once the number of samples goes above > 10^7.
You must be using a machine from several years ago, or are using an algorithm that requires lots of extra memory. 50*10^6 double- precision samples is under 400MB of memory. 10^8 double-precision numbers requires under a GB of memory, which should be fine on any recent computer. FFTW on my computer (which is what is used by Matlab) can compute an FFT of 10^8 double-precision real samples in under half a second on an Intel Core 2. You can calculate the autocorrelation of 10^8 samples with an FFT of only twice the length, so you should be fine. The advantage of using C instead of Matlab here would be that C gives you much better control over memory allocation; you might just be being careless and not realizing how many copies of your array you are making in Matlab. Regards, Steven G. Johnson
Reply by Rune Allnor April 17, 20092009-04-17
On 17 Apr, 10:21, SG <s.gesem...@gmail.com> wrote:
> On 17 Apr., 09:45, Rune Allnor <all...@tele.ntnu.no> wrote: > > > On 16 Apr, 23:27, SG <s.gesem...@gmail.com> wrote: > > > This statement is far too general. > > > Unless you have stripped the operations down to the core > > BLAS/LAPACK level, Vladimir is right. > > That's what I've been saying. Matlab comes with lots of built-in > functions that operate on vectors and matrices (find, min, max, > mldivide, +, -, *, .....) which is why I said "different programming > style" and "not easy to compare". No point in arguming here.
At least C++ (I don't know 'bout C) have a few libraries that provide data structures like trees, lists, queues etc, that allow for far more efficient algorithms than matlab can, which is locked into the 'straightjacket' posed by the ND array data structure. Here is a teaser I wrote some time ago on comp.soft-sys.matlab, about these very issues: http://groups.google.no/group/comp.soft-sys.matlab/msg/0b6492cb37b5d7fd?hl=no The main discussion was about some semantic detail about how to make matlab run fast. As you can see from the quoted profiler report, my 'naive' matlab alternative beat the 'elegant' matlab alternative by almost a factor 2. Then I did the same job with C++, the only difference being that I exploited the available data structures. This naive C++ code, with the algorith tailored to the problem but not optimized by me in any other way[*], ran almost 10x faster than the fastest matlab alternative. As I said, unless you work with core BLAS/LAPACK types of problems, leave matlab out if run-time is a concern.
> > There is the interpreter > > overhead as well as the call-by-value function call overhead, > > which take a lot more time than people realize. > > [...] > > And you have to throw every good programming practice > > down the drains. One aspect of matlab is that variables > > are passed between functions by copy, not by reference, > > which means that local copies of variables are found > > all over the place. > > Pass-by-value semantics doesn't exclude the possibility of "copy-on- > write" optimizations. It's perfectly possible to write > > &#4294967295; &#4294967295;x = somefunc(x); > > and let somefunc modify its local "copy" without any actual copying > being involved. Now, whether Matlab does this behind the scenes is > beyond me.
I am pretty sure it does. While not as bad as worst-case, it still represents a significant overhead compared to C or C++.
> Anyhow, I wasn't saying Matlab is better than anything else in every > situation! Obviously it isn't. And it might not be the best tool for > the original poster except for testing some quickly-written prototype > algorithms on "toy problems".
That's exactly what matlab was designed to do. Despite appearances, that's all matlab is useful for.
> But Vladimir's response was too general > and too one-sided which deserved some correcting, IMHO.
Vladimir made two very concise statements: 1) "Don't use MATLAB for serious computations." 2) "A plain C++ program runs ~10 times faster." The first I fully support, since matlab is a beast both what 'serious' programming/development and run-time resource demands are concerned. If you try and make anything beyond prototypes or toy programs, you will have to work your butt off in return for less than mediocre performance. As demonstrated in the post quoted above. Vladimir's second statementis factually correct for just about anything but core BLAS/LAPACK applications. Rune [*] I compiled the C++ code as a MEX routine to be called from within matlab. Whatever optimizing compiler directives were invoked from the MEX compilation setupwe, were standardized and not modified by me.
Reply by Rune Allnor April 17, 20092009-04-17
On 17 Apr, 09:54, itakatz <itak...@gmail.com> wrote:
> > 4) &#4294967295;Compute the DFT of the zero-padded frame > > 5) &#4294967295;Compute the periodogram as the squared magnitude > > &#4294967295; &#4294967295; of the spectrum coefficients > > 6) &#4294967295;Compute the IDFT of the periodogram > > > Rune > > A little bit off-topic: > > Since multiplication in frequency domain is the same as (circular) > convolving in the time domain,
Convolution by DFT is circular, which is why you need to zero-pad before computing the DFT. The convolution of two N-length sequences gives as result a sequence with at most 2N-1 non-zero elements. Zero-padding to at least 2N-1 total length allocates the space needed to avoid wrap-around effects. Depending on what implementation of the FFT you have available, you might zero-pad further, to the next power of 2.
> shouldn't one of the segments be time- > reversed before zero-padding and DFTing?
Time-reversing in time domain is not sufficient to handle wrap-around effects from circular convolution, since you don't change the length of the input sequence. With pen and paper your trick might work, but in the computer you need to allocate memory and find a map from the 'true' index space k = -N,...,N to the RAM address space, BaseAddress + {0,1,...2N}. Essentially, you need to do some book-keeping to find out if the rxx[0] coefficient turns out in the index 0 position of the buffer, or in the index N position.
> in that way the auto- > correlation is re-written as a convolution, and then it is ok to > compute it in the freq domain. Am I right?
In time domain, reversing the direction of a sequence x[n] is written as x[-n]. It's a standard excercise in most DSP books to show that if x[n] <-> X[k] is a DFT pair, then x[-n] <-> X'[k] where X' is the complex conjugate of X, is also a DFT pair. Well, this is correct at least when x[n] is real-valued, it might not be for complex-valued x[n]. Anyway, convolving x[n] by x[-n] in time domain, which is a common estimator for the autocorrelation function (save some scaling coefficients), can alternatively be computed in frequency domain as a product of the DFTs: Sxx[k] = X[k]X'[k] = |X[k]|^2. Hence my algorithm works, once the relevant scaling coefficients have been accounted for. Rune
Reply by SG April 17, 20092009-04-17
On 17 Apr., 09:45, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 16 Apr, 23:27, SG <s.gesem...@gmail.com> wrote: > > This statement is far too general. > > Unless you have stripped the operations down to the core > BLAS/LAPACK level, Vladimir is right.
That's what I've been saying. Matlab comes with lots of built-in functions that operate on vectors and matrices (find, min, max, mldivide, +, -, *, .....) which is why I said "different programming style" and "not easy to compare". No point in arguming here.
> There is the interpreter > overhead as well as the call-by-value function call overhead, > which take a lot more time than people realize. > [...] > And you have to throw every good programming practice > down the drains. One aspect of matlab is that variables > are passed between functions by copy, not by reference, > which means that local copies of variables are found > all over the place.
Pass-by-value semantics doesn't exclude the possibility of "copy-on- write" optimizations. It's perfectly possible to write x = somefunc(x); and let somefunc modify its local "copy" without any actual copying being involved. Now, whether Matlab does this behind the scenes is beyond me. Anyhow, I wasn't saying Matlab is better than anything else in every situation! Obviously it isn't. And it might not be the best tool for the original poster except for testing some quickly-written prototype algorithms on "toy problems". But Vladimir's response was too general and too one-sided which deserved some correcting, IMHO. Cheers! SG
Reply by itakatz April 17, 20092009-04-17
> 4) &#4294967295;Compute the DFT of the zero-padded frame > 5) &#4294967295;Compute the periodogram as the squared magnitude > &#4294967295; &#4294967295; of the spectrum coefficients > 6) &#4294967295;Compute the IDFT of the periodogram
> > Rune
A little bit off-topic: Since multiplication in frequency domain is the same as (circular) convolving in the time domain, shouldn't one of the segments be time- reversed before zero-padding and DFTing? in that way the auto- correlation is re-written as a convolution, and then it is ok to compute it in the freq domain. Am I right?
Reply by Rune Allnor April 17, 20092009-04-17
On 16 Apr, 23:27, SG <s.gesem...@gmail.com> wrote:
> On 16 Apr., 22:52, Vladimir Vassilevsky <antispam_bo...@hotmail.com> > wrote: > > > sasuke wrote: > > > Hi > > > > I am using MATLAB to generate close to 10^8 samples > > > First mistake. Don't use MATLAB for serious computations. A plain C++ > > program runs ~10 times faster. > > This statement is far too general.
Unless you have stripped the operations down to the core BLAS/LAPACK level, Vladimir is right. There is the interpreter overhead as well as the call-by-value function call overhead, which take a lot more time than people realize. More below.
> It's not easy to compare the > performances on this level. The programming style is different. If you > want your Matlab code to run fast you usually exploit all the built-in > functions (which map to ATLAS, FFTW, ....) and operate on a matrix/ > vector level as opposed to touching every scalar by hand.
And you have to throw every good programming practice down the drains. One aspect of matlab is that variables are passed between functions by copy, not by reference, which means that local copies of variables are found all over the place. This first of all consumes memory but also consumes lots of time copying variables back and forth, where one would use pointers or references in C or C++. The solution is to avoid using functions, thus forefeiting testing habits, code reusability and source code managment.
> But this may > also lead to higher memory requirements than what would be necessary > in case of C or C++.
Nope. It *might* not spend more memory. It *certainly* makes the program blow up and grab all the memory it can get. Rune
Reply by SG April 16, 20092009-04-16
On 16 Apr., 22:52, Vladimir Vassilevsky <antispam_bo...@hotmail.com>
wrote:
> sasuke wrote: > > Hi > > > I am using MATLAB to generate close to 10^8 samples > > First mistake. Don't use MATLAB for serious computations. A plain C++ > program runs ~10 times faster.
This statement is far too general. It's not easy to compare the performances on this level. The programming style is different. If you want your Matlab code to run fast you usually exploit all the built-in functions (which map to ATLAS, FFTW, ....) and operate on a matrix/ vector level as opposed to touching every scalar by hand. But this may also lead to higher memory requirements than what would be necessary in case of C or C++.
> > I am using MATLAB to do this. > > Don't use MATLAB. It is a toy for stupidents and studiots.
Man, improve your social skills. It may be of no great use to you. That doesn't mean it's not useful. You're basically insulting every Matlab user for no good reason. Matlab/Octave is nice for quickly testing things / prototyping / plotting 2D- and 3D- graphs. You can quickly import data and manipulate it. It has a nice collection of tools for all sorts of things and an interactive shell ... etc etc etc. C++ is nice. But it certainly isn't the best language for every job and you know it. SG