Power Spectral Density

In layman’s terms, the Power Spectral Density (PSD) shows the strength of the variations of the property that is measured as a function of frequency. In other words, it shows at which frequencies variations are strong and at which frequencies variations are weak 1 . The PSD will be used to normalize the power at all frequencies so that excess power at any frequency is more obvious which is a technique known as data whitening 2 .

The PSD is calculated by splitting up the signal into overlapping segments and scan through each segment to calculate individual periodogram. All the periodograms are then averaged, reducing the variance of the individual power measurements.

The length of the overall power spectral density can be determined from the maximum frequency and the segment’s length. Indeed, the maximum frequency determines the frequency range covered by the PSD while the inverse of the segment’s length gives the frequency separation between samples as it corresponds to the minimum frequency of a detectable signal. For instance, if a segment’s length is equal to 60 seconds and the maximum frequency in the filter bank is 256 Hz, the total length of the PSD is equal to 256 times 60 where 1/60 is the minimum frequency that can be detected. This can be programmed as follows:

# Minimum frequency of detectable signal in a segment
delta_f = 1. / psd_segment_length
# Calculate PSD length counting the zero frequency element
fd_len = fmax / delta_f + 1

In case an impulse response test is being made, a flat (white) PSD is created using the pycbc.types.FrequencySeries module as follows:

# Produce flat data
flat_data = numpy.ones(fd_len) * 2. / fd_len
# Create PSD frequency series
fd_psd = types.FrequencySeries(flat_data, 1./psd_segment_length, ts_data.start_time)

If no impulse response is considered, we then use the Welch’s method to perform the power spectral density estimate using the pycbc.psd.welch module:

# Create overall PSD using Welch's method
fd_psd = psd.welch(data,avg_method=psd_estimation,seg_len=seg_len,seg_stride=seg_stride)

What this will do is to compute the discrete Fourier transform for each PSD segment to produce invidual periodograms, and then compute the squared magnitude of the result. The individual periodograms are then averaged using the user-defined average method, psd_estimation , and return the frequency series, fd_psd , which will store the power measurement for each frequency bin.


A more detailed description of the PSD construction using the Welch’s method is provided in Section VI of Allan et al. (2012) .

If one wish, it is possible to display the power measurements, frequencies and frequency separation between consecutive samples as follows:

print 'Power measurements of the first 10 frequency bins:',
print fd_psd[:10]
print 'Frequencies of the first 10 bins:',
print fd_psd.sample_frequencies[:10]
print 'Frequency separation between bins:',
print fd_psd.delta_f

Finally, and in order for the fd_psd variable to be properly readable by the other LIGO tools while creating the filter bank, we need to convert it as follows:

# We need this for the SWIG functions
lal_psd = fd_psd.lal()