Search All of the Math Forum:

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

Notice: We are no longer accepting new posts, but the forums will continue to be readable.

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

 Messages: [ Previous | Next ]
 Jim Schwendeman Posts: 7 Registered: 6/21/13
Re: random time-series displacement data from PSD and back again
Posted: Jun 24, 2013 4:12 PM

> > 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:

% 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?

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