Time-Frequency Analysis

class Spectrogram(Array2D):
    def __init__(self, data, times=None, frequencies=None, t0=None, f0=None,
                 dt=None, df=None, unit=None, **kwargs):
        """
        Create a new Spectrogram.
        
        Parameters:
        - data: 2D array, spectrogram data (time x frequency)
        - times: array-like, time coordinates
        - frequencies: array-like, frequency coordinates  
        - t0: float, start time
        - f0: float, start frequency
        - dt: float, time resolution
        - df: float, frequency resolution
        - unit: str or Unit, data unit
        """
    
    @classmethod
    def read(cls, source, format=None, **kwargs):
        """
        Read Spectrogram from file.
        
        Parameters:
        - source: str, file path
        - format: str, file format
        
        Returns:
        Spectrogram object
        """
    
    def plot(self, **kwargs):
        """
        Create a plot of the spectrogram.
        
        Returns:
        Plot object with time-frequency visualization
        """
    
    def q_transform(self, **kwargs):
        """
        Calculate Q-transform of the spectrogram data.
        
        Returns:
        Q-transform Spectrogram
        """
    
    def ratio(self, method='median', **kwargs):
        """
        Calculate ratio against background.
        
        Parameters:
        - method: str, background estimation method
        
        Returns:
        Normalized Spectrogram
        """
    
    def crop(self, start=None, end=None, flow=None, fhigh=None, copy=True):
        """
        Crop to specific time-frequency region.
        
        Parameters:
        - start: float, start time
        - end: float, end time  
        - flow: float, low frequency
        - fhigh: float, high frequency
        - copy: bool, return copy or view
        
        Returns:
        Cropped Spectrogram
        """
    
    def zpk(self, zeros, poles, gain, **kwargs):
        """
        Apply frequency-domain filter.
        
        Parameters:
        - zeros: array-like, filter zeros
        - poles: array-like, filter poles
        - gain: float, filter gain
        
        Returns:
        Filtered Spectrogram
        """
    
    def percentile(self, percentile, **kwargs):
        """
        Calculate percentile across time axis.
        
        Parameters:  
        - percentile: float, percentile value (0-100)
        
        Returns:
        FrequencySeries with percentile values
        """
    
    def variance(self, **kwargs):
        """
        Calculate variance across time axis.
        
        Returns:
        FrequencySeries with variance values
        """

class SpectrogramList(list):
    def __init__(self, *args, **kwargs):
        """
        List of Spectrogram objects.
        """
    
    def join(self, gap='raise', **kwargs):
        """
        Join spectrograms into single Spectrogram.
        
        Parameters:
        - gap: str, gap handling ('raise', 'ignore', 'pad')
        
        Returns:
        Single joined Spectrogram
        """
    
    def plot(self, **kwargs):
        """
        Plot all spectrograms.
        
        Returns:
        Plot object
        """

from gwpy.signal.qtransform import q_scan

def q_scan(data, **kwargs):
    """
    Perform Q-transform scan for transient detection.
    
    Parameters:
    - data: TimeSeries, input time series data
    - qrange: tuple, Q-factor range
    - frange: tuple, frequency range  
    - mismatch: float, maximum mismatch
    - outseg: Segment, output time segment
    
    Returns:
    Spectrogram with Q-transform data
    """

from gwpy.timeseries import TimeSeries

# Read gravitational-wave strain data
strain = TimeSeries.fetch_open_data('H1', 1126259446, 1126259478,
                                   sample_rate=4096)

# Create spectrogram
spec = strain.spectrogram(stride=1,      # 1 second time bins
                         fftlength=4,    # 4 second FFTs  
                         overlap=2,      # 2 second overlap
                         window='hann')  # Hann window

# Plot spectrogram
plot = spec.plot(norm='log', vmin=1e-24, vmax=1e-19)
plot.set_xlabel('Time [s]')
plot.set_ylabel('Frequency [Hz]')
plot.set_title('H1 Strain Spectrogram')
plot.colorbar(label='Amplitude Spectral Density [strain/√Hz]')
plot.show()

from gwpy.signal import q_scan

# Focus on GW150914 time
gw_time = 1126259462.4
strain_zoom = strain.crop(gw_time-2, gw_time+2)

# Perform Q-transform scan
qspec = q_scan(strain_zoom, 
               qrange=(4, 64),           # Q-factor range
               frange=(20, 512),         # Frequency range
               mismatch=0.2,             # Maximum mismatch
               outseg=(gw_time-1, gw_time+1))

# Plot Q-transform
plot = qspec.plot(figsize=(12, 6))
plot.set_xlabel('Time [s]')
plot.set_ylabel('Frequency [Hz]') 
plot.set_title('Q-transform of GW150914')
plot.colorbar(label='Normalized energy')

# Mark GW time
plot.axvline(gw_time, color='red', linestyle='--', 
            label='GW150914')
plot.legend()
plot.show()

# Calculate background-normalized spectrogram
spec_ratio = spec.ratio(method='median')

# Plot ratio spectrogram to highlight transients
plot = spec_ratio.plot(vmin=0.1, vmax=10, norm='log')
plot.set_xlabel('Time [s]')
plot.set_ylabel('Frequency [Hz]')
plot.set_title('Normalized Spectrogram (Glitch Detection)')
plot.colorbar(label='Ratio to median')
plot.show()

# Find loud pixels for glitch identification
loud_times = spec_ratio.times[spec_ratio.max(axis=0) > 5]
print(f"Potential glitch times: {loud_times}")

# Calculate spectral percentiles over time
median_spectrum = spec.percentile(50)  # Median spectrum
p90_spectrum = spec.percentile(90)     # 90th percentile
p10_spectrum = spec.percentile(10)     # 10th percentile

# Plot spectral evolution
plot = median_spectrum.plot(label='Median', color='blue')
plot.add_frequencyseries(p90_spectrum, label='90th percentile', 
                        color='red', alpha=0.7)
plot.add_frequencyseries(p10_spectrum, label='10th percentile',
                        color='green', alpha=0.7)

plot.set_xlabel('Frequency [Hz]')
plot.set_ylabel('Amplitude Spectral Density [strain/√Hz]')
plot.set_yscale('log')
plot.legend()
plot.show()

from gwpy.timeseries import TimeSeriesDict

# Read data from multiple detectors
data = TimeSeriesDict.fetch_open_data(['H1', 'L1'], 
                                     start=1126259446, 
                                     end=1126259478)

# Create spectrograms for each detector
specs = {}
for ifo, ts in data.items():
    specs[ifo] = ts.spectrogram(stride=1, fftlength=4, overlap=2)

# Plot comparison
fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

for i, (ifo, spec) in enumerate(specs.items()):
    im = axes[i].imshow(spec.T, aspect='auto', origin='lower',
                       extent=[spec.xspan[0], spec.xspan[1],
                              spec.yspan[0], spec.yspan[1]],
                       norm='log', vmin=1e-24, vmax=1e-19)
    axes[i].set_ylabel('Frequency [Hz]')
    axes[i].set_title(f'{ifo} Spectrogram')
    
axes[1].set_xlabel('Time [s]')
plt.tight_layout()
plt.show()

# Calculate coherence spectrogram between detectors
h1_strain = TimeSeries.fetch_open_data('H1', start, end)
l1_strain = TimeSeries.fetch_open_data('L1', start, end)

# Coherence spectrogram  
coh_spec = h1_strain.coherence_spectrogram(l1_strain,
                                          stride=30,
                                          fftlength=4,
                                          overlap=2)

# Plot coherence
plot = coh_spec.plot(vmin=0, vmax=1)
plot.set_xlabel('Time [s]')
plot.set_ylabel('Frequency [Hz]')
plot.set_title('H1-L1 Coherence Spectrogram')
plot.colorbar(label='Coherence')
plot.show()

from gwpy.signal.qtransform import QTiling, QPlane

# Custom Q-transform with specific tiling
qtiling = QTiling(qrange=(4, 64), frange=(20, 512), mismatch=0.2)

# Calculate Q-transform manually for more control
qgram = strain_zoom.q_transform(qrange=(8, 32),
                               frange=(30, 400), 
                               gps=gw_time,
                               search=0.5)

# Plot with custom styling
plot = qgram.plot(figsize=(10, 6))
plot.set_ylim(30, 400)
plot.set_xlabel('Time from %.1f [s]' % gw_time)
plot.set_ylabel('Frequency [Hz]')
plot.set_title('Custom Q-transform')
plot.colorbar(label='Normalized energy')
plot.show()