FFT Spectrum Analyzer (python, written by Claude Code)

Claude Code wrote this python project from the task below. It installed and ran on a clean Ubuntu 22.04 machine as written.

Task given to the agent:

TASK: FFT Spectrum Analyzer

Write a program that performs Fast Fourier Transform (FFT) analysis on time-domain signal data, identifying dominant frequencies, computing power spectral density, and supporting windowing functions.

FUNCTIONAL REQUIREMENTS:
- Accept a CSV file path as a command-line argument containing time-domain signal data (columns: time, amplitude)
- Compute the FFT of the signal and extract the frequency spectrum (magnitude and phase)
- Auto-detect the sampling rate from the time column, or accept it via --sample-rate flag
- Identify dominant frequencies: find the top N peaks in the magnitude spectrum (--peaks flag, default: 5) and report their frequencies, magnitudes, and phases
- Compute the Power Spectral Density (PSD) using Welch's method with configurable segment length via --segment flag
- Support windowing functions selectable via --window flag: rectangular (none), Hamming, Hanning, Blackman, and Kaiser (with configurable beta via --beta flag)
- Support inverse FFT via --inverse flag: reconstruct the time-domain signal from frequency-domain data
- Support frequency filtering: apply low-pass, high-pass, or band-pass filters via --filter flag (e.g., --filter low:1000 for 1kHz low-pass) and output the filtered signal
- Export the frequency spectrum data as CSV via --export flag
- Print analysis summary to console: sampling rate, number of samples, frequency resolution, dominant frequencies with magnitudes, and total signal power
- Save the full analysis as JSON with --output flag (default: fft_analysis.json)
- If no input is given, generate a sample signal composed of three sine waves at known frequencies (100Hz, 250Hz, 800Hz) with added white noise, sampled at 4000Hz for 1 second, analyze it, and show that the FFT correctly identifies the three component frequencies
- Handle errors: non-uniform sampling, insufficient data points, and signals with DC offset

Create a complete Python project for a clean Ubuntu 22.04 machine with only Python 3.10+ installed. Include:
- Source code
- requirements.txt with all dependencies (direct and transitive) pinned to exact versions
- README.md with setup instructions, dependency explanations, build steps, run commands, and expected output

[]

"""
FFT Spectrum Analyzer
Performs FFT analysis on signals: frequency spectrum, PSD, windowing functions,
and frequency filtering.
Dependencies: numpy 1.26.4, scipy 1.12.0, matplotlib 3.8.2
"""

import numpy as np
from scipy import signal as scipy_signal
from scipy.fft import fft, fftfreq, ifft
import matplotlib.pyplot as plt


class FFTSpectrumAnalyzer:
    """Comprehensive FFT spectrum analyzer for signal processing."""

    WINDOW_FUNCTIONS = {
        "hann": np.hanning,
        "hamming": np.hamming,
        "blackman": np.blackman,
        "bartlett": np.bartlett,
        "kaiser": lambda n: np.kaiser(n, beta=14),
    }

    def __init__(self, sample_rate=1000.0):
        """Initialize the FFT Spectrum Analyzer.

        Args:
            sample_rate: Sampling rate in Hz.
        """
        self.sample_rate = sample_rate
        self.signal_data = None
        self.time_axis = None

    def generate_test_signal(self, duration, frequencies, amplitudes, noise_level=0.0):
        """Generate a composite test signal with optional noise.

        Args:
            duration: Duration of the signal in seconds.
            frequencies: List of frequency components in Hz.
            amplitudes: List of amplitude values for each frequency.
            noise_level: Standard deviation of Gaussian noise to add.

        Returns:
            Tuple of (time_axis, signal_data).
        """
        num_samples = int(self.sample_rate * duration)
        self.time_axis = np.linspace(0, duration, num_samples, endpoint=False)
        self.signal_data = np.zeros(num_samples)

        for freq, amp in zip(frequencies, amplitudes):
            self.signal_data += amp * np.sin(2 * np.pi * freq * self.time_axis)

        if noise_level > 0:
            self.signal_data += np.random.normal(0, noise_level, num_samples)

        return self.time_axis, self.signal_data

    def compute_fft(self, signal_data=None, window="hann"):
        """Compute the FFT of the signal with optional windowing.

        Args:
            signal_data: Input signal array. Uses stored signal if None.
            window: Window function name to apply before FFT.

        Returns:
            Dictionary with frequency axis, magnitude spectrum, and phase spectrum.
        """
        if signal_data is not None:
            data = np.array(signal_data, dtype=float)
        elif self.signal_data is not None:
            data = self.signal_data.copy()
        else:
            raise ValueError("No signal data available. Generate or provide signal first.")

        n = len(data)

        if window and window in self.WINDOW_FUNCTIONS:
            win = self.WINDOW_FUNCTIONS[window](n)
            data = data * win

        fft_result = fft(data)
        frequencies = fftfreq(n, 1.0 / self.sample_rate)

        positive_mask = frequencies >= 0
        freqs = frequencies[positive_mask]
        magnitude = 2.0 / n * np.abs(fft_result[positive_mask])
        magnitude[0] /= 2.0  # DC component
        phase = np.angle(fft_result[positive_mask])

        return {
            "frequencies": freqs,
            "magnitude": magnitude,
            "phase": phase,
            "fft_raw": fft_result,
        }

    def compute_psd(self, signal_data=None, method="welch", nperseg=256):
        """Compute Power Spectral Density of the signal.

        Args:
            signal_data: Input signal array. Uses stored signal if None.
            method: PSD estimation method ('welch' or 'periodogram').
            nperseg: Number of samples per segment for Welch method.

        Returns:
            Dictionary with frequencies and PSD values.
        """
        data = signal_data if signal_data is not None else self.signal_data
        if data is None:
            raise ValueError("No signal data available.")

        if method == "welch":
            freqs, psd = scipy_signal.welch(
                data, fs=self.sample_rate, nperseg=min(nperseg, len(data))
            )
        elif method == "periodogram":
            freqs, psd = scipy_signal.periodogram(data, fs=self.sample_rate)
        else:
            raise ValueError(f"Unknown PSD method: {method}")

        return {
            "frequencies": freqs,
            "psd": psd,
            "psd_db": 10 * np.log10(psd + 1e-12),
        }

    def apply_window(self, signal_data=None, window_type="hann"):
        """Apply a window function to the signal.

        Args:
            signal_data: Input signal. Uses stored signal if None.
            window_type: Type of window function to apply.

        Returns:
            Dictionary with windowed signal and window function values.
        """
        data = signal_data if signal_data is not None else self.signal_data
        if data is None:
            raise ValueError("No signal data available.")

        n = len(data)
        if window_type in self.WINDOW_FUNCTIONS:
            window = self.WINDOW_FUNCTIONS[window_type](n)
        else:
            raise ValueError(
                f"Unknown window type: {window_type}. "
                f"Available: {list(self.WINDOW_FUNCTIONS.keys())}"
            )

        windowed = data * window
        return {
            "windowed_signal": windowed,
            "window": window,
            "window_type": window_type,
        }

    def frequency_filter(self, signal_data=None, filter_type="lowpass",
                         cutoff=100.0, order=5, cutoff_high=None):
        """Apply a frequency domain filter to the signal.

        Args:
            signal_data: Input signal. Uses stored signal if None.
            filter_type: Type of filter ('lowpass', 'highpass', 'bandpass', 'bandstop').
            cutoff: Cutoff frequency in Hz (low cutoff for bandpass/bandstop).
            order: Filter order for Butterworth design.
            cutoff_high: High cutoff frequency for bandpass/bandstop filters.

        Returns:
            Dictionary with filtered signal and filter characteristics.
        """
        data = signal_data if signal_data is not None else self.signal_data
        if data is None:
            raise ValueError("No signal data available.")

        nyquist = self.sample_rate / 2.0

        if filter_type in ("bandpass", "bandstop"):
            if cutoff_high is None:
                raise ValueError(f"{filter_type} filter requires cutoff_high parameter.")
            wn = [cutoff / nyquist, cutoff_high / nyquist]
        else:
            wn = cutoff / nyquist

        b, a = scipy_signal.butter(order, wn, btype=filter_type)
        filtered = scipy_signal.filtfilt(b, a, data)

        w, h = scipy_signal.freqz(b, a, worN=2048, fs=self.sample_rate)

        return {
            "filtered_signal": filtered,
            "filter_response_freqs": w,
            "filter_response_magnitude": np.abs(h),
            "filter_response_db": 20 * np.log10(np.abs(h) + 1e-12),
            "filter_type": filter_type,
            "cutoff": cutoff,
            "order": order,
        }

    def plot_spectrum(self, fft_result, title="Frequency Spectrum", save_path=None):
        """Plot the frequency spectrum.

        Args:
            fft_result: Result dictionary from compute_fft().
            title: Plot title.
            save_path: Path to save the figure. Shows interactively if None.
        """
        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

        ax1.plot(fft_result["frequencies"], fft_result["magnitude"], color="steelblue")
        ax1.set_xlabel("Frequency (Hz)")
        ax1.set_ylabel("Magnitude")
        ax1.set_title(f"{title} - Magnitude Spectrum")
        ax1.grid(True, alpha=0.3)

        ax2.plot(fft_result["frequencies"], np.degrees(fft_result["phase"]),
                 color="coral", alpha=0.7)
        ax2.set_xlabel("Frequency (Hz)")
        ax2.set_ylabel("Phase (degrees)")
        ax2.set_title(f"{title} - Phase Spectrum")
        ax2.grid(True, alpha=0.3)

        plt.tight_layout()
        if save_path:
            plt.savefig(save_path, dpi=150, bbox_inches="tight")
            print(f"Spectrum plot saved to: {save_path}")
        else:
            plt.show()
        plt.close()

    def plot_psd(self, psd_result, title="Power Spectral Density", save_path=None):
        """Plot the Power Spectral Density.

        Args:
            psd_result: Result dictionary from compute_psd().
            title: Plot title.
            save_path: Path to save the figure.
        """
        fig, ax = plt.subplots(figsize=(12, 5))
        ax.semilogy(psd_result["frequencies"], psd_result["psd"], color="darkgreen")
        ax.set_xlabel("Frequency (Hz)")
        ax.set_ylabel("PSD (V^2/Hz)")
        ax.set_title(title)
        ax.grid(True, alpha=0.3)
        plt.tight_layout()

        if save_path:
            plt.savefig(save_path, dpi=150, bbox_inches="tight")
            print(f"PSD plot saved to: {save_path}")
        else:
            plt.show()
        plt.close()

    def plot_comparison(self, original, filtered, filter_result, save_path=None):
        """Plot comparison of original vs filtered signal with filter response.

        Args:
            original: Original signal data.
            filtered: Filtered signal data.
            filter_result: Result dictionary from frequency_filter().
            save_path: Path to save the figure.
        """
        fig, axes = plt.subplots(3, 1, figsize=(12, 10))
        time = np.arange(len(original)) / self.sample_rate

        axes[0].plot(time, original, alpha=0.7, label="Original")
        axes[0].plot(time, filtered, alpha=0.9, label="Filtered")
        axes[0].set_xlabel("Time (s)")
        axes[0].set_ylabel("Amplitude")
        axes[0].set_title("Time Domain Comparison")
        axes[0].legend()
        axes[0].grid(True, alpha=0.3)

        orig_fft = self.compute_fft(original, window=None)
        filt_fft = self.compute_fft(filtered, window=None)
        axes[1].plot(orig_fft["frequencies"], orig_fft["magnitude"],
                     alpha=0.7, label="Original")
        axes[1].plot(filt_fft["frequencies"], filt_fft["magnitude"],
                     alpha=0.9, label="Filtered")
        axes[1].set_xlabel("Frequency (Hz)")
        axes[1].set_ylabel("Magnitude")
        axes[1].set_title("Frequency Domain Comparison")
        axes[1].legend()
        axes[1].grid(True, alpha=0.3)

        axes[2].plot(filter_result["filter_response_freqs"],
                     filter_result["filter_response_db"], color="red")
        axes[2].set_xlabel("Frequency (Hz)")
        axes[2].set_ylabel("Gain (dB)")
        axes[2].set_title(
            f"Filter Response ({filter_result['filter_type']}, "
            f"cutoff={filter_result['cutoff']} Hz, order={filter_result['order']})"
        )
        axes[2].grid(True, alpha=0.3)
        axes[2].set_ylim(-80, 5)

        plt.tight_layout()
        if save_path:
            plt.savefig(save_path, dpi=150, bbox_inches="tight")
            print(f"Comparison plot saved to: {save_path}")
        else:
            plt.show()
        plt.close()

    def plot_windows(self, n=256, save_path=None):
        """Plot all available window functions and their frequency responses.

        Args:
            n: Window length.
            save_path: Path to save the figure.
        """
        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))

        for name, func in self.WINDOW_FUNCTIONS.items():
            win = func(n)
            ax1.plot(win, label=name)
            fft_win = fft(win, 2048)
            magnitude_db = 20 * np.log10(np.abs(fft_win[:1024]) + 1e-12)
            magnitude_db -= magnitude_db.max()
            freqs = np.linspace(0, 0.5, 1024)
            ax2.plot(freqs, magnitude_db, label=name)

        ax1.set_xlabel("Sample")
        ax1.set_ylabel("Amplitude")
        ax1.set_title("Window Functions (Time Domain)")
        ax1.legend()
        ax1.grid(True, alpha=0.3)

        ax2.set_xlabel("Normalized Frequency")
        ax2.set_ylabel("Magnitude (dB)")
        ax2.set_title("Window Functions (Frequency Domain)")
        ax2.legend()
        ax2.grid(True, alpha=0.3)
        ax2.set_ylim(-120, 5)

        plt.tight_layout()
        if save_path:
            plt.savefig(save_path, dpi=150, bbox_inches="tight")
            print(f"Window plot saved to: {save_path}")
        else:
            plt.show()
        plt.close()


def main():
    """Demonstrate the FFT Spectrum Analyzer capabilities."""
    print("=" * 60)
    print("FFT Spectrum Analyzer - Demo")
    print("=" * 60)

    analyzer = FFTSpectrumAnalyzer(sample_rate=1000.0)

    # Generate a test signal: 50 Hz + 120 Hz + 200 Hz with noise
    frequencies = [50, 120, 200]
    amplitudes = [1.0, 0.5, 0.3]
    t, sig = analyzer.generate_test_signal(
        duration=1.0, frequencies=frequencies,
        amplitudes=amplitudes, noise_level=0.2
    )
    print(f"\nGenerated signal: {len(sig)} samples at {analyzer.sample_rate} Hz")
    print(f"  Components: {list(zip(frequencies, amplitudes))}")
    print(f"  Duration: {t[-1]:.3f} s")

    # Compute FFT with different windows
    print("\n--- FFT Analysis with Different Windows ---")
    for window_name in ["hann", "hamming", "blackman", "bartlett"]:
        result = analyzer.compute_fft(window=window_name)
        peak_indices = np.argsort(result["magnitude"])[-3:][::-1]
        peak_freqs = result["frequencies"][peak_indices]
        peak_mags = result["magnitude"][peak_indices]
        print(f"  Window: {window_name:10s} | Top peaks: ", end="")
        for f, m in zip(peak_freqs, peak_mags):
            print(f"{f:.1f} Hz ({m:.4f})", end="  ")
        print()

    # Power Spectral Density
    print("\n--- Power Spectral Density ---")
    psd_welch = analyzer.compute_psd(method="welch", nperseg=256)
    psd_per = analyzer.compute_psd(method="periodogram")
    print(f"  Welch PSD:       {len(psd_welch['frequencies'])} frequency bins")
    print(f"  Periodogram PSD: {len(psd_per['frequencies'])} frequency bins")
    print(f"  Welch peak freq: {psd_welch['frequencies'][np.argmax(psd_welch['psd'])]:.1f} Hz")
    print(f"  Periodogram peak: {psd_per['frequencies'][np.argmax(psd_per['psd'])]:.1f} Hz")

    # Windowing functions comparison
    print("\n--- Window Function Properties ---")
    for name, func in analyzer.WINDOW_FUNCTIONS.items():
        win = func(256)
        print(f"  {name:10s}: sum={np.sum(win):.2f}, "
              f"max={np.max(win):.4f}, "
              f"coherent_gain={np.mean(win):.4f}")

    # Frequency filtering
    print("\n--- Frequency Filtering ---")

    # Low-pass filter (remove 120 Hz and 200 Hz components)
    lp_result = analyzer.frequency_filter(filter_type="lowpass", cutoff=80.0, order=5)
    lp_fft = analyzer.compute_fft(lp_result["filtered_signal"], window=None)
    peak_idx = np.argmax(lp_fft["magnitude"])
    print(f"  Low-pass (80 Hz):  dominant freq = {lp_fft['frequencies'][peak_idx]:.1f} Hz")

    # High-pass filter (remove 50 Hz component)
    hp_result = analyzer.frequency_filter(filter_type="highpass", cutoff=90.0, order=5)
    hp_fft = analyzer.compute_fft(hp_result["filtered_signal"], window=None)
    peak_idx = np.argmax(hp_fft["magnitude"])
    print(f"  High-pass (90 Hz): dominant freq = {hp_fft['frequencies'][peak_idx]:.1f} Hz")

    # Band-pass filter (isolate 120 Hz component)
    bp_result = analyzer.frequency_filter(
        filter_type="bandpass", cutoff=100.0, cutoff_high=140.0, order=4
    )
    bp_fft = analyzer.compute_fft(bp_result["filtered_signal"], window=None)
    peak_idx = np.argmax(bp_fft["magnitude"])
    print(f"  Band-pass (100-140 Hz): dominant freq = {bp_fft['frequencies'][peak_idx]:.1f} Hz")

    # Band-stop filter (notch out 120 Hz component)
    bs_result = analyzer.frequency_filter(
        filter_type="bandstop", cutoff=110.0, cutoff_high=130.0, order=4
    )
    bs_fft = analyzer.compute_fft(bs_result["filtered_signal"], window=None)
    peak_idx = np.argmax(bs_fft["magnitude"])
    print(f"  Band-stop (110-130 Hz): dominant freq = {bs_fft['frequencies'][peak_idx]:.1f} Hz")

    # Summary
    print("\n" + "=" * 60)
    print("Analysis complete.")
    print(f"  Signal length:      {len(sig)} samples")
    print(f"  Sample rate:        {analyzer.sample_rate} Hz")
    print(f"  Frequency range:    0 - {analyzer.sample_rate / 2:.0f} Hz")
    print(f"  Frequency resolution: {analyzer.sample_rate / len(sig):.2f} Hz")
    print("=" * 60)


if __name__ == "__main__":
    main()

# FFT Spectrum Analyzer (Python - Trial 1)

Performs FFT analysis on signals: frequency spectrum, PSD, windowing functions, and frequency filtering.

## Dependencies

- **numpy** (1.26.4) - Numerical array operations and signal generation
- **scipy** (1.12.0) - FFT computation, signal processing, and filter design
- **matplotlib** (3.8.2) - Plotting frequency spectra and filter responses

## Setup

```bash
pip install -r requirements.txt
```

## Usage

```bash
python fft_analyzer.py
```

## Features

- Composite signal generation with configurable frequencies and amplitudes
- FFT computation with multiple window functions (Hann, Hamming, Blackman, Bartlett, Kaiser)
- Power Spectral Density estimation via Welch and periodogram methods
- Frequency filtering (low-pass, high-pass, band-pass, band-stop) using Butterworth design
- Visualization of spectra, PSD, filter responses, and window function comparisons

numpy==1.26.4
scipy==1.12.0
matplotlib==3.8.2