random time-series displacement data from PSD and back again

[Jim Schwendeman]

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

[TideMan]

Post by Jim Schwendeman
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...
% 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%
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.

[Jim Schwendeman]

Post by TideMan

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.

Hi Tideman,

Thanks for the response. I found the examples as suggested, and was able to duplicate fairly well, but still have a couple of questions. First, current code posted below:

%% Description:
% This generates random time-series displacement/velocity data from a PSD
% using the method described here:

% http://www.mathworks.com/matlabcentral/newsreader/view_thread/42194

% 1. Take the sqrt of your PSD. This gives the amplitude spectrum.
% 2. Generate uniform noise for the length of the spectrum using rand
% and multiply by 2pi. This is a set of random phases.
% 3. Make a set of complex noise using Y=amp.*exp(i*phi)
% 4. Replicate complex conj Y=[Y;flipud(conj(Y))];
% 5. Give it zero mean Y=[0;Y]
% 6. Inverse transform y=ifft(Y)
% 7. Scale to get the desired the std devn

function [t,x,xdot, Fs] = CreateTimeSeries(N, FLAG_plot)

%% Create original PSD

N = 4000
Fs = N
F_Max = 2000
gRMS = 2.5

freq = linspace(0,F_Max,N)';
PSD = ones(size(freq));
cumulative = cumtrapz(freq,PSD);

%Scale PSD to desired gRMS
PSD = (gRMS^2/max(cumulative)).*PSD;

if(FLAG_plot)
figure(11)
subplot(2,1,1)
semilogx(freq,PSD,'LineWidth',2)
hold on
end

%% Create ASD
ASD = sqrt(PSD);

%% Create random phases
rand_phase = 2*pi*rand(size(ASD));

%% Make complex noise
Y =ASD.*exp(i*rand_phase);

%% Replicate Complex Conj
Y=[Y;flipud(conj(Y))];

%% Give it 0 Mean:
Y=[0;Y];

%% IFFT
Fs = N
dt = 1/N

%Length doubled when we replicated complex conjugate above
t = 0:dt:2;
y=ifft(Y);

%Parceval's theorem works as follows:
% Sum(0:N-1) x(i)^2 = (1/N)SUM(F(i)^2) = Average g^2/Hz
% SO... sqrt(avg g^2/Hz)*sqrt(Hz) gets us the area = gRMS
gRMS_time_series = sqrt(sum(abs(y).^2))*sqrt(F_Max)

if(FLAG_plot)
figure(11)
subplot(2,1,2)
plot(t,y)
xlabel('time (s)')
ylabel('displacement (m)')
title(['gRMS_{time series} = ',num2str(gRMS_time_series)])
end

disp(['Average g^2/Hz = ',num2str(sum(abs(y).^2))])

%% CONVERT BACK TO FREQUENCY DOMAIN FOR GOOD FEELINGS
x = y;

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);
freq_new = 0:(2*F_Max)/length(x):F_Max;

psdx=psdx*(length(psdx))^2;

gRMS_recreate = sqrt(max(cumtrapz(freq_new,psdx)))

if(FLAG_plot)
figure(11)
subplot(2,1,1)
semilogx(freq_new,psdx,'r');
grid on;
title(['gRMS_{original} = ',num2str(gRMS_original),', gRMS_{recreate} = ',num2str(gRMS_recreate)])
legend('PSD_{orig}','PSD_{recreate}')
xlabel('Frequency (Hz)');
ylabel('g^2/Hz');
end

%Where I'm confused now is the relative scaling when going between the ifft and fft. When I converted from the time domain back to the frequency domain, I'm just taking half the data, multiplying by 2, and that gets me to the correct g^2/Hz amplitude. You'll note that I multiply and divide the psdx value by N^2, which doesn't make sense. I thought going between we'd always scale by sqrt(N).

Thoughts?

[Jim Schwendeman]

Post by TideMan

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.

Hi Tideman,

Thanks for the response. I found the examples as suggested, and was able to duplicate fairly well, but still have a couple of questions. First, current code posted below:

%% Description:
% This generates random time-series displacement/velocity data from a PSD
% using the method described here:

% http://www.mathworks.com/matlabcentral/newsreader/view_thread/42194

% 1. Take the sqrt of your PSD. This gives the amplitude spectrum.
% 2. Generate uniform noise for the length of the spectrum using rand
% and multiply by 2pi. This is a set of random phases.
% 3. Make a set of complex noise using Y=amp.*exp(i*phi)
% 4. Replicate complex conj Y=[Y;flipud(conj(Y))];
% 5. Give it zero mean Y=[0;Y]
% 6. Inverse transform y=ifft(Y)
% 7. Scale to get the desired the std devn

function [t,x,xdot, Fs] = CreateTimeSeries(N, FLAG_plot)

%% Create original PSD

N = 4000
Fs = N
F_Max = 2000
gRMS = 2.5

freq = linspace(0,F_Max,N)';
PSD = ones(size(freq));
cumulative = cumtrapz(freq,PSD);

%Scale PSD to desired gRMS
PSD = (gRMS^2/max(cumulative)).*PSD;

if(FLAG_plot)
figure(11)
subplot(2,1,1)
semilogx(freq,PSD,'LineWidth',2)
hold on
end

%% Create ASD
ASD = sqrt(PSD);

%% Create random phases
rand_phase = 2*pi*rand(size(ASD));

%% Make complex noise
Y =ASD.*exp(i*rand_phase);

%% Replicate Complex Conj
Y=[Y;flipud(conj(Y))];

%% Give it 0 Mean:
Y=[0;Y];

%% IFFT
Fs = N
dt = 1/N

%Length doubled when we replicated complex conjugate above
t = 0:dt:2;
y=ifft(Y);

%Parceval's theorem works as follows:
% Sum(0:N-1) x(i)^2 = (1/N)SUM(F(i)^2) = Average g^2/Hz
% SO... sqrt(avg g^2/Hz)*sqrt(Hz) gets us the area = gRMS
gRMS_time_series = sqrt(sum(abs(y).^2))*sqrt(F_Max)

if(FLAG_plot)
figure(11)
subplot(2,1,2)
plot(t,y)
xlabel('time (s)')
ylabel('displacement (m)')
title(['gRMS_{time series} = ',num2str(gRMS_time_series)])
end

disp(['Average g^2/Hz = ',num2str(sum(abs(y).^2))])

%% CONVERT BACK TO FREQUENCY DOMAIN FOR GOOD FEELINGS
x = y;

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);
freq_new = 0:(2*F_Max)/length(x):F_Max;

psdx=psdx*(length(psdx))^2;

gRMS_recreate = sqrt(max(cumtrapz(freq_new,psdx)))

if(FLAG_plot)
figure(11)
subplot(2,1,1)
semilogx(freq_new,psdx,'r');
grid on;
title(['gRMS_{original} = ',num2str(gRMS_original),', gRMS_{recreate} = ',num2str(gRMS_recreate)])
legend('PSD_{orig}','PSD_{recreate}')
xlabel('Frequency (Hz)');
ylabel('g^2/Hz');
end

%Where I'm confused now is the relative scaling when going between the ifft and fft. When I converted from the time domain back to the frequency domain, I'm just taking half the data, multiplying by 2, and that gets me to the correct g^2/Hz amplitude. You'll note that I multiply and divide the psdx value by N^2, which doesn't make sense. I thought going between we'd always scale by sqrt(N).

Thoughts?

[Jim Schwendeman]

TideMan,

Maybe the better question is this:

When I create my time-series data, how do I know that it's scaled correctly? When I replicated the complex conjugate prior to the ifft, I doubled the number of points, so I guess I'm not entirely sure.

I believe I am correct, because if I do the sum(abs(y).^2), that is equal to the mean g^2/Hz from the original PSD, which leads me to believe I've satisfied Parseval's theorem.

Thoughts?

[TideMan]

Post by Jim Schwendeman
TideMan,
When I create my time-series data, how do I know that it's scaled correctly? When I replicated the complex conjugate prior to the ifft, I doubled the number of points, so I guess I'm not entirely sure.
I believe I am correct, because if I do the sum(abs(y).^2), that is equal to the mean g^2/Hz from the original PSD, which leads me to believe I've satisfied Parseval's theorem.
Thoughts?

I can't follow your code. You seem to define PSD as ones, then use it as if it had data. I gave up.............

This is Parseval's Law:
sum(PSD)*df must equal var(y)
The area under the spectrum must be the variance of the time signal.
In your case, I think df is 1.