Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Topic: random time-series displacement data from PSD and back again
Replies: 5   Last Post: Jun 24, 2013 9:01 PM

 Messages: [ Previous | Next ]
 Derek Goring Posts: 3,922 Registered: 12/7/04
Re: random time-series displacement data from PSD and back again
Posted: Jun 21, 2013 5:48 PM

On Saturday, June 22, 2013 8:10:26 AM UTC+12, Jim Schwendeman wrote:
> Howdy Forum,
>
>
>
> I have a project where I have a vibration PSD profile (G^2/Hz). I'd like to use this to create random time-series displacement data to feed into a dynamic model to evaluate performance in the presence of vibe.
>
>
>
> The problem occurs when I convert from G^2/Hz to displacement/time, and then back to the identical (or very near it) G^2/Hz profile. I create random time-series data, but when I convert back to the frequency domain, it doesn't match up. The re-created PSD has content up to twice the frequency domain I'm initially putting in, and the RMS values don't match up. This sanity check needs to work before I can confidently put the time-series data into the model, so I'm a bit hung up...
>
>
>
> Code below:
>
> % Create noise (centered 1), pump it through PSD, and then calc gRMS
>
> noise_range = 0.02
>
> noise = noise_range*randn(1,N)+(1-noise_range/2);
>
> PSD_noise = PSD.*noise;
>
>
>
> %Calculate gRMS for original PSD and noise-injected PSD
>
> gRMS = sqrt(max(cumtrapz(PSD0(:,1),PSD0(:,2))));
>
> gRMS_noise = sqrt(max(cumtrapz(f,PSD_noise)));
>
> %These match up, the gRMS of my noisy PSD is nearly identical to my input PSD
>
>
>
> % Create time and displacement
>
> t = linspace(0, N*dt, N);
>
> x = (1/sqrt(N))*real(ifft(sqrt(PSD_noise)));
>
> x = (x*9.81)./sqrt(std(x));
>
> disp('Parceval"s Theorem says this result should = gRMS')
>
> sum(abs(x).^2)
>
> %This answer is similar to the gRMS of my noisy/input PSD, but now off by ~10%
>
>
>
> %Now I convert time-series displacement back to frequency domain and should get the same PSD back out:
>
> N = length(x);
>
> xdft = fft(x);
>
> xdft = xdft(1:N/2+1);
>
> psdx = (1/(Fs*N)).*abs(xdft).^2;
>
> psdx(2:end-1) = 2*psdx(2:end-1);
>
> psdx = psdx*sqrt(N)*9.81;
>
> freq = 0:Fs/length(x):Fs/2;
>
>
>
> When I end up plotting freq vs psdx, I get frequency content over double of what was put in, and my gRMS is orders of magnitude below what was expected.
>
>
>
> Any thoughts?
>
> Jim

If PSD is the traditional power spectral density, then you have not prepared the Fourier transform correctly before ifft.
To see how to do this correctly, search "phase randomisation" (note spelling) on this newsgroup.

Date Subject Author
6/21/13 Jim Schwendeman
6/21/13 Derek Goring
6/24/13 Jim Schwendeman
6/24/13 Jim Schwendeman
6/24/13 Jim Schwendeman
6/24/13 Derek Goring