Creating a time series signal from a known PSD

Question

alex bradley 2023-5-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1961044-creating-a-time-series-signal-from-a-known-psd

评论： Star Strider 2023-5-11

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

alex bradley 2023-5-11

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 2023-5-11

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 2023-5-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1961044-creating-a-time-series-signal-from-a-known-psd#answer_1232794

在 MATLAB Online 中打开

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)