Autocorrelation matrix of a long sequence

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

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.

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].

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