Forums

Finding psd and autocorrelation for a random number matrix 50x50

Started by akhil1449 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