PSD curve into time series data

78 次查看（过去 30 天）
Chandramouli Jambunathan
评论： 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更早的评论
Chandramouli Jambunathan
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
dpb 2022-7-8
That link is to a FEX (FileExchange) function; it gives examples of use there...

请先登录，再进行评论。

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

回答（2 个）

Chunru
Chunru 2022-7-8
Ran in:
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])
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.
  2 个评论
Chunru
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)

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by