Frequency Domain Analysis

class FrequencySeries(Series):
    def __init__(self, data, frequencies=None, f0=None, df=None, 
                 unit=None, name=None, **kwargs):
        """
        Create a new FrequencySeries.
        
        Parameters:
        - data: array-like, frequency series data values
        - frequencies: array-like, frequency values for each sample
        - f0: float, starting frequency (Hz)
        - df: float, frequency spacing (Hz)
        - unit: str or Unit, physical unit of the data
        - name: str, name for this series
        """
    
    @classmethod
    def read(cls, source, format=None, **kwargs):
        """
        Read FrequencySeries from file.
        
        Parameters:
        - source: str, path to data file
        - format: str, file format ('txt', 'hdf5', 'ligolw')
        
        Returns:
        FrequencySeries object
        """
    
    def plot(self, **kwargs):
        """
        Create a plot of the frequency series.
        
        Returns:
        Plot object
        """
    
    def zpk(self, zeros, poles, gain, **kwargs):
        """
        Filter using zero-pole-gain representation.
        
        Parameters:
        - zeros: array-like, filter zeros
        - poles: array-like, filter poles  
        - gain: float, filter gain
        
        Returns:
        Filtered FrequencySeries
        """
    
    def crop(self, start=None, end=None, copy=True):
        """
        Crop to specific frequency range.
        
        Parameters:
        - start: float, start frequency
        - end: float, end frequency
        - copy: bool, return copy or view
        
        Returns:
        Cropped FrequencySeries
        """
    
    def interpolate(self, df, **kwargs):
        """
        Interpolate to new frequency spacing.
        
        Parameters:
        - df: float, new frequency resolution
        
        Returns:
        Interpolated FrequencySeries
        """
    
    def inject(self, other, **kwargs):
        """
        Inject another frequency series.
        
        Parameters:  
        - other: FrequencySeries, series to inject
        
        Returns:
        Combined FrequencySeries
        """
    
    def write(self, target, format=None, **kwargs):
        """
        Write FrequencySeries to file.
        
        Parameters:
        - target: str, output file path
        - format: str, output format
        """

class SpectralVariance(Array2D):
    def __init__(self, data, bins=None, **kwargs):
        """
        2D histogram for spectral variance analysis.
        
        Parameters:
        - data: 2D array, spectral variance data
        - bins: array-like, histogram bin edges
        """
    
    @classmethod
    def from_spectrogram(cls, spectrogram, **kwargs):
        """
        Calculate spectral variance from spectrogram.
        
        Parameters:
        - spectrogram: Spectrogram, input data
        
        Returns:
        SpectralVariance object
        """
    
    def plot(self, **kwargs):
        """
        Plot the spectral variance.
        
        Returns:
        Plot object
        """

from gwpy.timeseries import TimeSeries
from gwpy.frequencyseries import FrequencySeries

# Calculate PSD from time series
strain = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478)
psd = strain.psd(fftlength=4, overlap=2)

# Plot PSD
plot = psd.plot()
plot.set_xlabel('Frequency [Hz]')
plot.set_ylabel('Power Spectral Density [strain^2/Hz]')
plot.set_yscale('log')
plot.show()

# Read PSD from file
ref_psd = FrequencySeries.read('reference_psd.txt')

# Compare PSDs
ratio = psd / ref_psd
ratio.plot()

# Calculate ASD (square root of PSD)
asd = strain.asd(fftlength=4)

# Plot with proper scaling
plot = asd.plot()
plot.set_xlabel('Frequency [Hz]')
plot.set_ylabel('Amplitude Spectral Density [strain/√Hz]')
plot.set_yscale('log')
plot.set_xlim(10, 2000)
plot.show()

from scipy.signal import butter

# Design filter in zpk format
z, p, k = butter(8, [30, 400], btype='bandpass', fs=4096, output='zpk')

# Apply filter to frequency series data
filtered_psd = psd.zpk(z, p, k)

# Crop to analysis band
analysis_psd = psd.crop(20, 2000)

# Read calibration transfer function
cal_response = FrequencySeries.read('calibration_response.txt', 
                                   format='txt')

# Apply calibration correction
calibrated_psd = raw_psd * cal_response

# Plot magnitude and phase
fig, axes = plt.subplots(2, 1, sharex=True)

# Magnitude
axes[0].loglog(cal_response.frequencies, abs(cal_response))
axes[0].set_ylabel('Magnitude')

# Phase  
axes[1].semilogx(cal_response.frequencies, 
                 np.angle(cal_response, deg=True))
axes[1].set_ylabel('Phase [degrees]')
axes[1].set_xlabel('Frequency [Hz]')

from gwpy.frequencyseries import SpectralVariance

# Calculate spectrogram first
spec = strain.spectrogram(stride=60, fftlength=4)

# Create spectral variance histogram
variance = SpectralVariance.from_spectrogram(spec, 
                                           bins=50, 
                                           low=1e-24, 
                                           high=1e-20)

# Plot spectral variance
plot = variance.plot()
plot.set_xlabel('Frequency [Hz]')
plot.set_ylabel('Amplitude Spectral Density [strain/√Hz]')
plot.show()

# Read multiple PSDs for noise budget
psds = {}
channels = ['H1:CAL-DELTAL_EXTERNAL_DQ',
           'H1:ASC-AS_A_RF45_I_ERR_DQ', 
           'H1:LSC-DARM_IN1_DQ']

for channel in channels:
    ts = TimeSeries.get(channel, start, end)
    psds[channel] = ts.psd(fftlength=8)

# Plot noise budget
plot = FrequencySeries.plot(figsize=(12, 6))
for name, psd in psds.items():
    plot.add_frequencyseries(psd, label=name)

plot.set_xlabel('Frequency [Hz]')
plot.set_ylabel('Power Spectral Density') 
plot.set_yscale('log')
plot.legend()
plot.show()

# Calculate cross-spectral density between channels
witness = TimeSeries.get('H1:ASC-AS_A_RF45_I_ERR_DQ', start, end)
strain = TimeSeries.get('H1:DCS-CALIB_STRAIN_C02', start, end)

# Cross-spectral density
csd = strain.csd(witness, fftlength=4, overlap=2)

# Coherence calculation
coherence = strain.coherence(witness, fftlength=4, overlap=2)

# Plot coherence
plot = coherence.plot()
plot.set_xlabel('Frequency [Hz]')
plot.set_ylabel('Coherence')
plot.set_ylim(0, 1)
plot.show()