# Finding psd and autocorrelation for a random number matrix 50x50

Started by August 12, 2015
```HI, I am a mechanical engineer, to generate surface roughness, I am
following an algorithm. I need to find autocorrelation and psd of the
input random matrix.

v=randn(50,50)
tempvar = xcorr2(v,v);
[size_1 size_2] = size(tempvar);
size_1 = (size_1 + 1)/2 - 1;
R = (toeplitz(tempvar(size_1+1 : 2*size_1+1)))./size_1;
psd=abs(fftshift(fft(R)));

By theory, the psd must be a constant value, but I get a range of values
varying from 0.02 to 3.7

Kindly advise me if there is any problem in this code?

---------------------------------------
Posted through http://www.DSPRelated.com
```
```On Wed, 12 Aug 2015 08:54:23 -0500, akhil1449 wrote:

> HI, I am a mechanical engineer, to generate surface roughness, I am
> following an algorithm. I need to find autocorrelation and psd of the
> input random matrix.
>
> v=randn(50,50)
> tempvar = xcorr2(v,v);
> [size_1 size_2] = size(tempvar);
> size_1 = (size_1 + 1)/2 - 1;
> R = (toeplitz(tempvar(size_1+1 : 2*size_1+1)))./size_1;
> psd=abs(fftshift(fft(R)));
>
>
> By theory, the psd must be a constant value, but I get a range of values
> varying from 0.02 to 3.7
>
> Kindly advise me if there is any problem in this code?

I'm not sure what all you're doing, yet I think I know what your problem
is.

Power spectral density is the expected value of the signal power at any
given time and frequency.  It is NOT the ACTUAL signal power for some
finite-time, real-world sample.

I don't know what quirks the 2D surface throws at you, but to get the PSD
of a 1D sample (the time series of a signal is what I'm used to), you
need to compute FFTs for a large number of same-sized segments of the
signal, calculate the power, and then average all of those power
measurements together to get an estimate of the PSD.

Off the top of my head I don't know how many segments you need to average
together to get a good result for PSD -- but I do know that the number
depends on how you define "good": the less deviation between your
measured number and the underlying actual number that you can stand, the
more segments you need to average.

--

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
```
```Tim Wescott  <seemywebsite@myfooter.really> wrote:

>On Wed, 12 Aug 2015 08:54:23 -0500, akhil1449 wrote:

>> HI, I am a mechanical engineer, to generate surface roughness, I am
>> following an algorithm. I need to find autocorrelation and psd of the
>> input random matrix.

>> v=randn(50,50)

[snip]

>> By theory, the psd must be a constant value, but I get a range of values
>> varying from 0.02 to 3.7
>>
>> Kindly advise me if there is any problem in this code?
>
>I'm not sure what all you're doing, yet I think I know what your problem
>is.
>
>Power spectral density is the expected value of the signal power at any
>given time and frequency.  It is NOT the ACTUAL signal power for some
>finite-time, real-world sample.
>
>I don't know what quirks the 2D surface throws at you, but to get the PSD
>of a 1D sample (the time series of a signal is what I'm used to), you
>need to compute FFTs for a large number of same-sized segments of the
>signal, calculate the power, and then average all of those power
>measurements together to get an estimate of the PSD.
>
>Off the top of my head I don't know how many segments you need to average
>together to get a good result for PSD -- but I do know that the number
>depends on how you define "good": the less deviation between your
>measured number and the underlying actual number that you can stand, the
>more segments you need to average.

I agree with this.  More generally, in calculating a PSD, the
autocorrelation function must be windowed, that is some sort
of ensemble average is being performed.  It is not just a straight
function over the entire dataset, and the size of the dataset
must be "larger" than the length of the window, otherwise, the
results will make little sense.

"Larger" could be a factor of 10 in many cases.

By not windowing, instead of seeing a flat spectrum typical
of white noise, you are seeing that spectrum convolved with
that of something basically rectangular.  So you expect some
slope-off from the center to the edges.

Offhand, applying a Hamming window to your xcorr2() output
might work.  You may need to start with bigger than 50 x 50.

Good luck

Steve
```
```On Thursday, August 13, 2015 at 1:54:27 AM UTC+12, akhil1449 wrote:
> HI, I am a mechanical engineer, to generate surface roughness, I am
> following an algorithm. I need to find autocorrelation and psd of the
> input random matrix.
>
> v=randn(50,50)
> tempvar = xcorr2(v,v);
> [size_1 size_2] = size(tempvar);
> size_1 = (size_1 + 1)/2 - 1;
> R = (toeplitz(tempvar(size_1+1 : 2*size_1+1)))./size_1;
> psd=abs(fftshift(fft(R)));
>
>
> By theory, the psd must be a constant value, but I get a range of values
> varying from 0.02 to 3.7
>
> Kindly advise me if there is any problem in this code?
>
>
>
>
>
> ---------------------------------------
> Posted through http://www.DSPRelated.com

assuming a rectangular window you will need more than just the one matrix but a collection of say 128 such matrices.
```
```<gyansorova@gmail.com> wrote:

>> v=randn(50,50)
>> tempvar = xcorr2(v,v);
>> [size_1 size_2] = size(tempvar);
>> size_1 = (size_1 + 1)/2 - 1;
>> R = (toeplitz(tempvar(size_1+1 : 2*size_1+1)))./size_1;
>> psd=abs(fftshift(fft(R)));

>> Kindly advise me if there is any problem in this code?

>assuming a rectangular window you will need more than just the one
>matrix but a collection of say 128 such matrices.

Depending on what you're trying to achieve, you may or may not
want to do this.  One can find the PSD of a single dataset,
or a collection of datasets, and either might be indicated
by the circumstances.

If you did want to collect 128 such matrices, in the above
code snippet, you would call xcorr2(v,v) 128 times on the
128 datasets, and average the results.

But this doesn't solve the problem mentioned before, that
the resulting autocorrelation function is never windowed.

Steve
```
```On Monday, August 17, 2015 at 5:05:54 AM UTC+12, Steve Pope wrote:
> <gyansorova@gmail.com> wrote:
>
> >> v=randn(50,50)
> >> tempvar = xcorr2(v,v);
> >> [size_1 size_2] = size(tempvar);
> >> size_1 = (size_1 + 1)/2 - 1;
> >> R = (toeplitz(tempvar(size_1+1 : 2*size_1+1)))./size_1;
> >> psd=abs(fftshift(fft(R)));
>
> >> Kindly advise me if there is any problem in this code?
>
> >assuming a rectangular window you will need more than just the one
> >matrix but a collection of say 128 such matrices.
>
> Depending on what you're trying to achieve, you may or may not
> want to do this.  One can find the PSD of a single dataset,
> or a collection of datasets, and either might be indicated
> by the circumstances.
>
> If you did want to collect 128 such matrices, in the above
> code snippet, you would call xcorr2(v,v) 128 times on the
> 128 datasets, and average the results.
>
> But this doesn't solve the problem mentioned before, that
> the resulting autocorrelation function is never windowed.
>
> Steve

Sure, but from a sampled data set you cannot get consistent results with an FFT without averaging.
```
```<gyansorova@gmail.com> wrote:

>On Monday, August 17, 2015 at 5:05:54 AM UTC+12, Steve Pope wrote:

>> Depending on what you're trying to achieve, you may or may not
>> want to do this.  One can find the PSD of a single dataset,
>> or a collection of datasets, and either might be indicated
>> by the circumstances.

>> If you did want to collect 128 such matrices, in the above
>> code snippet, you would call xcorr2(v,v) 128 times on the
>> 128 datasets, and average the results.

>> But this doesn't solve the problem mentioned before, that
>> the resulting autocorrelation function is never windowed.

>Sure, but from a sampled data set you cannot get consistent results with
>an FFT without averaging.

I am not sure what you're getting at.

The PSD is the cosine transform of the autocorrelation function.
The autocorellation function needs to come from somewhere, and
to be statistically valid needs to be formed from sufficient data
for one's purposes, but that sufficient data could just as easily be
a single long sequence instead of multiple shorter sequences.

So, to repeat, there is nothing inherent about forming a PSD
that requires any sort of ensemble averaging before forming
the autocorrelation.

As for FFT's, an FFT function can be used to compute some (not all)
cosine transforms, and that is its only role here.

Hope this makes some sense.

Steve

```