ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Thu, 02 Aug 2012 07:11:13 +0200Spectral density of a signalhttps://ask.sagemath.org/question/9148/spectral-density-of-a-signal/ Hello!
Could somebody please help me find a way to build the spectral density function for a given signal in Sage?
Thank you!Tue, 17 Jul 2012 09:48:47 +0200https://ask.sagemath.org/question/9148/spectral-density-of-a-signal/Answer by v_2e for <p>Hello!</p>
<p>Could somebody please help me find a way to build the spectral density function for a given signal in Sage?</p>
<p>Thank you!</p>
https://ask.sagemath.org/question/9148/spectral-density-of-a-signal/?answer=13886#post-id-13886A self-answer:
The **p**ower **s**pectral **d**ensity (PSD) [may be defined](https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density) as
$ S(\omega) = \lim \limits_{T \to +\infty} \frac{\left \vert F_T(\omega) \right \vert ^2}{T} $,
where $ F_T (\omega)$ is the Fourier transform defined as follows:
$ {F}_T(\omega) = \int \limits_0^T f(t) \exp(-i\omega t) ~ dt$
def PSD(time_series):
import scipy
signal_length = n(len(time_series)*(time_series[1][0]-time_series[0][0]))
signal_fft = scipy.fft(zip(*time_series)[1])
spectrum = []
for i in range(len(signal_fft)//2):
spectrum.append((i/signal_length,abs(signal_fft[i])**2/len(time_series)))
return spectrum
The accepted data set format is:
data = [(t1,y1),(t2,y2),...,(tn,yn)]
Calling
PSD(data)
for such signal will return the Power Spectral Density of a signal.
----------
Sometimes it is useful to apply some kind of [window function](https://en.wikipedia.org/wiki/Window_function) to a signal prior to calculating the PSD, since the sharp start and end of the data record may produce some spurious spectral components.
Here is an example of the popular [Hanning Window](https://en.wikipedia.org/wiki/Window_function#Hann_.28Hanning.29_window) application for the time series:
def hanning_window(time_series):
''' Applies Hanning window to the time series.
Accepted data format is a list of tuples [(x1,y1),(x2,y2),...]'''
series_length = n(len(time_series))
processed_signal=[]
for i in range(len(time_series)):
processed_signal.append((time_series[i][0], n(time_series[i][1] * \
(0.5*(1 - cos(2*pi*i/(series_length-1)))))))
return processed_signal
The result of its application looks like this:
![The result of Hanning window function influcence on the signal](/upfiles/1435950346516497.png)
One can simply call
PSD(hanning_window(data))
to get the power spectral density for a data set with Hanning window function applied.
You can compare the results of spectral density calculation for the initial time series and "windowed" time series:
![image description](/upfiles/1435950416656873.png)
Thu, 02 Aug 2012 07:11:13 +0200https://ask.sagemath.org/question/9148/spectral-density-of-a-signal/?answer=13886#post-id-13886