PSD curve into time series data
95 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (2)
Chunru
le 8 Juil 2022
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 commentaires
Chandramouli Jambunathan
le 11 Juil 2022
2 commentaires
Chunru
le 12 Juil 2022
Modifié(e) : Chunru
le 12 Juil 2022
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)
Voir également
Catégories
En savoir plus sur Fourier Analysis and Filtering dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!