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: phase randomisation, Zero Padding of IFFT, and Parseval
Replies: 2   Last Post: Jul 3, 2013 5:13 PM

 Messages: [ Previous | Next ]
 Jim Schwendeman Posts: 7 Registered: 6/21/13
phase randomisation, Zero Padding of IFFT, and Parseval
Posted: Jul 3, 2013 10:30 AM

Howdy forum,

My goal is to take a PSD (g^2/Hz), and generate random time series acceleration data of arbitrary length and sampling rate. As recommended in a different post, I am utilizing phase randomisation, which seems to be working as expected.

The process I'm using is as follows:
1.) Sample the PSD at the desired sampling rate from 0 to F_Max and convert to ASD
2.) zero pad ASD to length N = Total_time(say 5 sec)*Fs
3.) ifftshift
4.) phase randomisation
5.) At this point, I check that Parseval's is satisfied, meaning that the mean_squared_ASD == sum_squared_time_series_accel.
6.) I then run the FFT on the time series data to reconstruct the PSD, which also gives me confidence that the time-series data was scaled appropriately.

Questions that are giving me a bit (lot) of trouble:
1.) How do I know that the time series I've created (should be 0:5 secs), is correct? Parseval's doesn't take into account time duration, so while the amplitude comparison is correct, how do I know I haven't mistakenly extended or compressed time?

2.) For zero padding prior to the IFFT, does it matter if I zeropad the center, vs zeropad the outsides and then call IFFTSHIFT?

I've attached code that should run w/o any dependencies:

clc
clear all

%% Define Params
PSD0 = [0 1000 2000;
0.0035 0.0035 0.0023]';

Fs = 1000
T = 5
N = T*Fs

%% Calcs

F_Max = max(PSD0(:,1))
gRMS_original = sqrt(max(cumtrapz(PSD0(:,1),PSD0(:,2))))

freq = linspace(0,F_Max,Fs)';
PSD = interp1(PSD0(:,1),PSD0(:,2),freq);

figure(11)
subplot(2,1,1)
semilogx(freq,PSD,'LineWidth',2)
hold on
xlabel('Hz')
ylabel('PSD g^2/Hz')

% So NOW have our PSD, sampled at some known frequency, now zeropad
aSD = sqrt(PSD);

sum_squared_aSD = sum(abs(zp_aSD).^2)/(N-1)

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

%% Make complex noise
A=zp_aSD.*exp(i*rand_phase);

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

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

%% Have to ditch NaN
A(isnan(A)) = 1;
A(isinf(A)) = 1;

%% IFFT
dt = T/N

%Length doubled when we replicated complex conjugate above
t = 0:dt:T - dt;
a=ifft(A, N);

figure(11)
subplot(2,1,2)
plot(t,a)

disp(['Verify Parseval"s Theorem'])
sum_accel_tseries = sum(abs(a.^2))

x = a;
fft_x = fft(x,N);
xdft = abs(ifftshift(fft_x));

%Have to square to go from aSD to PSD
psdx = (xdft).^2;

freq_new = 0:(F_Max)/Fs:F_Max - (F_Max/Fs);

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

subplot(2,1,1)
semilogx(freq_new,psdx,'r');
title(['gRMS_{recreate} = ',num2str(gRMS_recreate)])
grid on;

error = sum_squared_aSD-sum_accel_tseries

Date Subject Author
7/3/13 Jim Schwendeman
7/3/13 Jim Schwendeman
7/3/13 Jim Schwendeman