Creating a time series signal from a known PSD

61 visualizaciones (últimos 30 días)
alex bradley
alex bradley el 10 de Mayo de 2023
Comentada: Star Strider el 11 de Mayo de 2023
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 comentarios
alex bradley
alex bradley el 11 de Mayo de 2023
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
Star Strider el 11 de Mayo de 2023
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.

Iniciar sesión para comentar.

Respuestas (1)

Chunru
Chunru el 11 de Mayo de 2023
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)
[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 fs/2])

Categorías

Más información sobre Signal Generation and Preprocessing en Help Center y File Exchange.

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by