Is there a code to get power spectrum corresponding to a velocity/acceleration time history?

I am working on generating seimic noise numerically. I have a velocity/acceleration time history. I want to generate power spectrum (power spectral density) compatible to the time history. Is there any code for that?

[Sam Chak]

Hi @Aniket Desai

Best if you provide the noise data and the equation for the power spectrum, so that we can investigate and analyze it.

Else, please check if the following article is helpful:

https://www.mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

[Aniket Desai]

I have attached the noise data. Left column is time (s) and right column is velocity (m/s).

Answer 1

[William Rose]

@Aniket Desai,

[edited: corrected spelling errors]

Yes it is possible. Estimate the spectrum of the original signal by taking the FFT or by other standard methods, such as those described in the link provided by @Sam Chak, or see the attached notes. Note that the Welch periodogram method, which is usually my first choice, would not be my first choice for a seismogram, because Welch's method chops up the signal into segments that are presumed to be similar random samples of the signal. But for a seismogram, they are not random samples of the signal, because one chunk will have the peak, and other chunks will have very little activity, and so on. Therefore I would just use the fft to estimate the amplitude spectrum.

Then create a frequency domain representation of your signal, with the ampltudes you got above, and random phases. But the phases for the "negative" frequencies must be the negative of the phases at the corresponding positive frequencies, otherwise the inverse FFT will be complex, and you want it to be real. Then take the inverse FFT of the frequency-domain version of the signal which you have cretaed. This is the output you desire. I willl provide an example of all this, in case it is not clear.

[William Rose]

I am attaching a synthetic seismogram, not very realistic, to illustrate the points I made above. It is a single column of 200 numbers. Let us assume that these were sampled at 25 Hz. Therefore the signal is 8 seconds long. Will continue later, gotta run.

[William Rose]

Here is an example of how to generate signals with the same power spectrum as the original signal. A simulated signal is attached, different than before: now it is 4 seconds sampled at 100 Hz: seismogramSynth.txt. The earlier signal had a different duration, sampling rate, and envelope. The script reads in the data from a real (or simulated, in this case) seismogram. Then it synthesizes three seismograms with the same power at every frequency, but with random phase at each frequency. The script includes explanatory comments.

Top left plot: Simulated seismogram. Bottom left: Amplitude spectrum of simulated seismogram, up to the Nyquist frequency.

Top right: Three synthetic seismograms. Bottom right: Amplitude spectra of simulated seismograms.

The amplitude spectra in the figure show that the synthetic seismograms (lower right plot) have the exact same power at each frequency as the signal which was read in from a file (lower left plot). However, the phases of the 3 synthetic seismograms are different, and in fact are random. The synthetic seismograms look wrong: they don't have an envelope that quickly rises and then gradually decays, like real earthquake seismograms. The relative phases of a real earthquake seismogram are NOT random - they are correlated in a way that produces the typical envelope. By randomizing the phases, we have destroyed the shape of the envelope. But the three synthetic signals do have what you requested: a power spectrum that matches the power spectrum of the original.

[Sam Chak]

@William Rose Good job! 👍

[William Rose]

@Sam Chak, thank you very much.

[William Rose]

@Aniket Desai,

I modified the script so that it reads in your file and uses the time vector in your file. I also modified the code so it works with an even or odd number of points - since your file has an odd number of points. As before, the script synthesizes three signals with the same power spectrum as the original signal, but with random phases. If you wish, adjust the script to save the signals to disk, or change the number of signals. Here is the graphical output, using your data file.