Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



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 timeseries 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); zerosPAD = (NFs)/2 zp_aSD = ifftshift([zeros(zerosPAD,1); aSD; zeros(zerosPAD,1)]);
sum_squared_aSD = sum(abs(zp_aSD).^2)/(N1)
%% 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)); xdft = xdft(zerosPAD + 1:length(xdft)zerosPAD);
%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_aSDsum_accel_tseries



