Creating a time series signal from a known PSD

I am trying to use a PSD I have generated (which represents a frequency response function) to filter (white noise) time series data. I usually use a software called ncode which uses what they call a "Custom Fourier Filter" to do this. However, for this application it would be advantageous to use matlab. I have cobbled together some code using guesswork and google however, any thoughts on how to go about this would be very welcome. I have attached a .mfile here. Thanks in advance for your help!

I should also mention I am not sure if fft or ifft is the right way round!

f = [0,0.0732422000000000,0.146484000000000,0.219727000000000,0.292969000000000,0.366211000000000,0.439453000000000,0.512695000000000,0.585938000000000,0.659180000000000,0.732422000000000,0.805664000000000,0.878906000000000,0.952148000000000,1.02539000000000,1.09863000000000,1.17188000000000,1.24512000000000,1.31836000000000,1.39160000000000,1.46484000000000,1.53809000000000,1.61133000000000,1.68457000000000,1.75781000000000,1.83105000000000,1.90430000000000,1.97754000000000,2.05078000000000,2.12402000000000,2.19727000000000,2.27051000000000,2.34375000000000,2.41699000000000,2.49023000000000,2.56348000000000,2.63672000000000,2.70996000000000,2.78320000000000,2.85645000000000,2.92969000000000,3.00293000000000,3.07617000000000,3.14941000000000,3.22266000000000,3.29590000000000,3.36914000000000,3.44238000000000,3.51563000000000,3.58887000000000,3.66211000000000,3.73535000000000,3.80859000000000,3.88184000000000,3.95508000000000,4.02832000000000,4.10156000000000,4.17480000000000,4.24805000000000,4.32129000000000,4.39453000000000,4.46777000000000,4.54102000000000,4.61426000000000,4.68750000000000,4.76074000000000,4.83398000000000,4.90723000000000,4.98047000000000]; % need f(1)=0; f(end)=fs/2

p = [0.685300000000000,0.569000000000000,0.391500000000000,0.960200000000000,0.715100000000000,0.710800000000000,0.541300000000000,0.656100000000000,0.669000000000000,0.768500000000000,0.794200000000000,0.789000000000000,0.756500000000000,0.986900000000000,0.921900000000000,0.796700000000000,0.890700000000000,1.04290000000000,1.07830000000000,1.06320000000000,1.06600000000000,0.973400000000000,0.665000000000000,0.694600000000000,0.759300000000000,0.798900000000000,0.932900000000000,0.913900000000000,0.965800000000000,0.785400000000000,0.838700000000000,0.861600000000000,0.881000000000000,0.994600000000000,0.974900000000000,0.851100000000000,0.865100000000000,0.789100000000000,0.889700000000000,1.10680000000000,1.03190000000000,0.995900000000000,1.28190000000000,1.17300000000000,1.18160000000000,1.37010000000000,1.45800000000000,1.08010000000000,1.16140000000000,1.05390000000000,1.05510000000000,1.20390000000000,1.06880000000000,0.936000000000000,0.978000000000000,1.14230000000000,1.27420000000000,0.822600000000000,1.38300000000000,1.74600000000000,2.92540000000000,1.30420000000000,0.632000000000000,0.997600000000000,1.72530000000000,1.12590000000000,0.783400000000000,1.25600000000000,1.11690000000000]; % power

fs = 10; % 2x fmax

% interpolate the spectrum

nfft = 4096;

f1 = linspace(0, fs/2, nfft)';

p1 = interp1(f, p, f1); % any suitable interp method

% design FIR filter which has desired response

hfir = fir2(nfft, f1/(fs/2), sqrt(p1)); % p1 is power

[ph, fh] =freqz(hfir, 1, 8192, fs);

hfir1 = hfir * sqrt(fs/2); % normalized freq -> freq

% generate time series

x = filter(hfir1, 1, randn(nfft*8, 1)); % generate data of length nfft*8

[px, fx] = pwelch(x, nfft, 3*nfft/4, nfft, fs); % estimate spectrum

plot(fh, 20*log10(abs(ph)), 'r-', f, 10*log10(p), 'bo', f1, 10*log10(p1), 'k:', fx, 10*log10(px), 'g-')

xlabel('f(Hz)'); ylabel('dB'); legend('Designed', 'Specified', 'Interp', 'Sig');

xlim([0 500])

[Star Strider]

A time series requires both the magnitude and phase in order to invert the Fourier transform. The PSD calculation destroys the imaginary part of the original complex Fourier transform, so inverting it to create a time series is not possible.

[alex bradley]

I should also mention I do have the complex data but, google mashing means I didn't include it in the code. Are there any useful examples of this?

[alex bradley]

Also I am not completely interested in recounstructing the exact time series but rather a time series composed of the same frequency content. I am not sure but, doesn't that mean I don't nessisarily need the imaginary part?

[Star Strider]

If you have the complex data, and if they are for the full, two-sided Fourier transform, just use ifft. (If you used fftshift, ise ifftshift first). If you have a one-sided Fourier transform, reconstruct the second half by flipping (flip) the conjugate (conj) of the original, and concantenating it to the end of the existing vector.

The sampling interval will be the same as the original data. If you only have the Nyquist frequency ‘Fn’, the sampling interval is:

Ts = 1/(2*Fn);

That should allow you to regererate the time vector.

Answer 1

[Chunru]

f = [0,0.0732422000000000,0.146484000000000,0.219727000000000,0.292969000000000,0.366211000000000,0.439453000000000,0.512695000000000,0.585938000000000,0.659180000000000,0.732422000000000,0.805664000000000,0.878906000000000,0.952148000000000,1.02539000000000,1.09863000000000,1.17188000000000,1.24512000000000,1.31836000000000,1.39160000000000,1.46484000000000,1.53809000000000,1.61133000000000,1.68457000000000,1.75781000000000,1.83105000000000,1.90430000000000,1.97754000000000,2.05078000000000,2.12402000000000,2.19727000000000,2.27051000000000,2.34375000000000,2.41699000000000,2.49023000000000,2.56348000000000,2.63672000000000,2.70996000000000,2.78320000000000,2.85645000000000,2.92969000000000,3.00293000000000,3.07617000000000,3.14941000000000,3.22266000000000,3.29590000000000,3.36914000000000,3.44238000000000,3.51563000000000,3.58887000000000,3.66211000000000,3.73535000000000,3.80859000000000,3.88184000000000,3.95508000000000,4.02832000000000,4.10156000000000,4.17480000000000,4.24805000000000,4.32129000000000,4.39453000000000,4.46777000000000,4.54102000000000,4.61426000000000,4.68750000000000,4.76074000000000,4.83398000000000,4.90723000000000,4.98047000000000];                 % need f(1)=0; f(end)=fs/2
p = [0.685300000000000,0.569000000000000,0.391500000000000,0.960200000000000,0.715100000000000,0.710800000000000,0.541300000000000,0.656100000000000,0.669000000000000,0.768500000000000,0.794200000000000,0.789000000000000,0.756500000000000,0.986900000000000,0.921900000000000,0.796700000000000,0.890700000000000,1.04290000000000,1.07830000000000,1.06320000000000,1.06600000000000,0.973400000000000,0.665000000000000,0.694600000000000,0.759300000000000,0.798900000000000,0.932900000000000,0.913900000000000,0.965800000000000,0.785400000000000,0.838700000000000,0.861600000000000,0.881000000000000,0.994600000000000,0.974900000000000,0.851100000000000,0.865100000000000,0.789100000000000,0.889700000000000,1.10680000000000,1.03190000000000,0.995900000000000,1.28190000000000,1.17300000000000,1.18160000000000,1.37010000000000,1.45800000000000,1.08010000000000,1.16140000000000,1.05390000000000,1.05510000000000,1.20390000000000,1.06880000000000,0.936000000000000,0.978000000000000,1.14230000000000,1.27420000000000,0.822600000000000,1.38300000000000,1.74600000000000,2.92540000000000,1.30420000000000,0.632000000000000,0.997600000000000,1.72530000000000,1.12590000000000,0.783400000000000,1.25600000000000,1.11690000000000];     % power
fs = 10;      % 2x fmax
% interpolate the spectrum
nfft = 4096;
f1 = linspace(0, fs/2, nfft)';

% need right interp/extrap method so that there is no nan in p1 
p1 = interp1(f, p, f1, 'linear', 'extrap')     % any suitable interp method

p1 = 4096×1

0.6853 0.6834 0.6814 0.6795 0.6775 0.6756 0.6737 0.6717 0.6698 0.6679

% any(isnan(p1))   % previously p1 contains nans
% design FIR filter which has desired response
%whos
hfir = fir2(nfft, f1/(fs/2), sqrt(p1))  % p1 is power

hfir = 1×4097

1.0e+00 * -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000

plot(hfir)

[ph, fh] =freqz(hfir, 1, 8192, fs);
hfir1 = hfir * sqrt(fs/2);        % normalized freq -> freq
% generate time series
x = filter(hfir1, 1, randn(nfft*8, 1));              % generate data of length nfft*8
[px, fx] = pwelch(x, nfft, 3*nfft/4, nfft, fs);        % estimate spectrum
plot(fh, 20*log10(abs(ph)), 'r-', f, 10*log10(p), 'bo', f1, 10*log10(p1), 'k:', fx, 10*log10(px), 'g-')
xlabel('f(Hz)'); ylabel('dB'); legend('Designed', 'Specified', 'Interp', 'Sig');
xlim([0 fs/2])