PSD curve into time series data

Question

Chandramouli Jambunathan 2022-7-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1754630-psd-curve-into-time-series-data

评论： Chandramouli Jambunathan 2022-7-14

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?

3 个评论
显示 1更早的评论隐藏 1更早的评论

Chandramouli Jambunathan 2022-7-7

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 2022-7-8

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Chunru 2022-7-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1754630-psd-curve-into-time-series-data#answer_1002685

在 MATLAB Online 中打开

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])

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Chandramouli Jambunathan 2022-7-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1754630-psd-curve-into-time-series-data#answer_1004595

@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.

2 个评论
显示无隐藏无

Chunru 2022-7-12

编辑：Chunru 2022-7-12

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 2022-7-14