PSD curve into time series data

Hello , I am trying to convert the PSD curve into time series data.

I came across this link https://uk.mathworks.com/matlabcentral/fileexchange/47342-timeseriesfrompsd-sxx-fs-t-plot_on, but its not clear what is the format of Sxx.

All I have is a 2X2 matrix of PSD curve

example:

Freq(Hz) 10 30 100 120 500 2000

PSD(G2/Hz) .01 .002 .002 .0004 .0004 .0004

Can somebody help me?

[dpb]

The FEX function expects you to provide it a 2-sided PSD; the DC in the first position and then the frequency compmonents to Nyquist followed by the reversed...

Something like

P2=[0 PSD fliplr(PSD)]/2;

where I've assumed a zero DC component.

With no more data than those, you're not likely to get much of a time realization, but it should match the PSD data if compute the PSD from the result which is all the info you've got.

[Chandramouli Jambunathan]

Hi @dpb,

Thanks for your quick response.

Sorry I didn't understand enough. What is FEX function? Is it an inbuilt function in matlab?currently I am using 2018a matlab version.

How do I give the below table to the function?

Freq(Hz) 10 30 100 120 500 2000

PSD(G2/Hz) .01 .002 .002 .0004 .0004 .0004

I want output time series sample rate at 60 Hz(60 samples per second) and total length of time as 5 min or more.

Thanks for your support.

[dpb]

That link is to a FEX (FileExchange) function; it gives examples of use there...

Answer 1

[Chunru]

f = [0 10 30 100 120 500 2000];                 % need f(1)=0; f(end)=fs/2
p = [0.01 .01 .002 .002 .0004 .0004 .0004];     % power

fs = 4000;      % 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])

Answer 2

[Chandramouli Jambunathan]

@Chunru Thanks a lot for the code and useful comments.

Is the signal 'ph' in frequency domain?

I am trying to convert 'ph' to time-domain using ifft (inverse fft), but initial results shows unrealistic acceleration values( Y- axis).

My objective is to plot acceleration(unit G) Vs time(unit seconds or minutes)

Now I am using another known PSD Vibration profile to double check the code

f = [10 28 40 100 200 500 2000]; in hertz

p = [0.02 0.02 0.04 0.04 0.08 0.08 0.00505]; in G^2/Hz

In time-domain, I am expecting the Y-axis(Acceleration in G) to be 7.94 GRMS(Root mean square).

Please let me know if you have any thoughts.

Note: 'interp1' resulted in NaN values, so I used 'pchip' interpolation method.

Thanks for your time and effort.

[Chunru]

The principle of generating the random noise with the specified spectra in above code is to filter the wihite noise (which has a psd of constant through out freq) via a spectum-shaping filter (hfir in the code). So we start with the spectrum specification and designs a FIR filte (it could be other type of filter, with different design method). In above, we use fir2 filter design method to get the filter coefficients hfir.

> Is the signal 'ph' in frequency domain?

'ph' is the filter's frequency response. It is in the frequency domain.

> I am trying to convert 'ph' to time-domain using ifft (inverse fft), but initial results shows unrealistic acceleration values( Y- axis).

'ph' is frequency response of the filter. IFFT of ph is the filter inpulse response hfir. You need to filter white noise through the filter hfir to get the realistic acceleration.

> My objective is to plot acceleration(unit G) Vs time(unit seconds or minutes)

The random time series data is x = filter(hfir1, 1, randn(nfft*8, 1))

> Now I am using another known PSD Vibration profile to double check the code

> f = [10 28 40 100 200 500 2000]; in hertz

> p = [0.02 0.02 0.04 0.04 0.08 0.08 0.00505]; in G^2/Hz

> In time-domain, I am expecting the Y-axis(Acceleration in G) to be 7.94 GRMS(Root mean square).

plot(x) or rms(x) to check if the random signal generated meet your requirement.

> Note: 'interp1' resulted in NaN values, so I used 'pchip' interpolation method.

If you use interp1, you must specify the spectrum value at f=0 and f=fs/2. (you have option of extrap, but it may not be robust)

[Chandramouli Jambunathan]

Thanks @Chunru.

I will look into this.