main.tex

\documentclass[11pt,a4paper]{article}

% Signal Processing and Fourier Analysis Template
% SEO Keywords: signal processing latex template, fourier analysis latex, digital signal processing latex

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{mathtools}
\usepackage{siunitx}
\usepackage{graphicx}
\usepackage{float}
\usepackage{hyperref}
\usepackage{cleveref}
\usepackage{booktabs}
\usepackage{pythontex}

\usepackage[style=numeric,backend=biber]{biblatex}
\addbibresource{references.bib}

\title{Signal Processing and Fourier Analysis:\\
Digital Methods and Spectral Techniques}
\author{CoCalc Scientific Templates}
\date{\today}

\begin{document}

\maketitle

\begin{abstract}
This comprehensive signal processing template demonstrates Fourier analysis, digital filter design, spectral estimation, and time-frequency analysis. Features include FFT algorithms, window functions, filter implementations, and wavelet transforms with applications to real-world signal analysis problems and professional visualizations.

\textbf{Keywords:} signal processing, Fourier analysis, digital filters, spectral estimation, FFT, wavelets, time-frequency analysis
\end{abstract}

\tableofcontents
\newpage

\section{Fourier Transform and Spectral Analysis}
\label{sec:fourier-analysis}

\begin{pycode}
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
from scipy import signal
from scipy.signal import windows
from scipy.fft import fft, fftfreq, fftshift

# Set random seed for reproducibility
np.random.seed(42)

def generate_test_signal():
    """Generate complex test signal for analysis"""

    # Time parameters
    fs = 1000  # Sampling frequency (Hz)
    T = 2.0    # Duration (seconds)
    t = np.arange(0, T, 1/fs)
    N = len(t)

    # Signal components
    f1, f2, f3 = 50, 120, 300  # Frequencies (Hz)

    # Construct signal: sinusoids + chirp + noise
    signal_clean = (np.sin(2*np.pi*f1*t) +
                   0.5*np.sin(2*np.pi*f2*t + np.pi/4) +
                   0.3*np.sin(2*np.pi*f3*t))

    # Add frequency-modulated chirp
    chirp = signal.chirp(t, f0=10, f1=200, t1=T, method='linear')

    # Add noise
    noise = 0.2 * np.random.randn(N)

    test_signal = signal_clean + 0.3*chirp + noise

    print(f"Signal Processing: Test Signal Generation")
    print(f"  Sampling frequency: {fs} Hz")
    print(f"  Duration: {T} s")
    print(f"  Samples: {N}")
    print(f"  Component frequencies: {f1}, {f2}, {f3} Hz")

    return t, test_signal, signal_clean, fs

def fourier_analysis(t, x, fs):
    """Perform comprehensive Fourier analysis"""

    N = len(x)

    # Compute FFT
    X = fft(x)
    freqs = fftfreq(N, 1/fs)

    # Magnitude and phase spectra
    magnitude = np.abs(X)
    phase = np.angle(X)

    # Power spectral density
    psd = magnitude**2 / (fs * N)

    # Only positive frequencies for display
    positive_freqs = freqs[:N//2]
    magnitude_positive = magnitude[:N//2]
    psd_positive = psd[:N//2]

    print(f"\nFourier Analysis Results:")
    print(f"  Frequency resolution: {freqs[1]:.2f} Hz")
    print(f"  Maximum frequency: {positive_freqs[-1]:.1f} Hz")

    # Find spectral peaks
    peaks, properties = signal.find_peaks(magnitude_positive,
                                        height=np.max(magnitude_positive)*0.1,
                                        distance=10)

    peak_freqs = positive_freqs[peaks]
    peak_magnitudes = magnitude_positive[peaks]

    print(f"  Detected peaks at frequencies:")
    for freq, mag in zip(peak_freqs, peak_magnitudes):
        print(f"    {freq:.1f} Hz (magnitude: {mag:.0f})")

    return freqs, X, magnitude, phase, psd, positive_freqs, magnitude_positive, psd_positive

# Generate and analyze test signal
t, test_sig, clean_sig, fs = generate_test_signal()
freqs, X, mag, phase, psd, pos_freqs, mag_pos, psd_pos = fourier_analysis(t, test_sig, fs)
\end{pycode}

\section{Digital Filter Design}
\label{sec:digital-filters}

\begin{pycode}
# Digital filter design and implementation
def design_filters(fs):
    """Design various digital filters"""

    nyquist = fs / 2

    # Low-pass Butterworth filter
    cutoff_lp = 100  # Hz
    order_lp = 4
    b_lp, a_lp = signal.butter(order_lp, cutoff_lp/nyquist, btype='low')

    # High-pass Butterworth filter
    cutoff_hp = 200  # Hz
    order_hp = 4
    b_hp, a_hp = signal.butter(order_hp, cutoff_hp/nyquist, btype='high')

    # Band-pass Chebyshev filter
    low_band, high_band = 80, 150  # Hz
    order_bp = 6
    b_bp, a_bp = signal.cheby1(order_bp, 1, [low_band/nyquist, high_band/nyquist], btype='band')

    # FIR filter using window method
    numtaps = 101
    cutoff_fir = 150  # Hz
    b_fir = signal.firwin(numtaps, cutoff_fir/nyquist, window='hamming')
    a_fir = [1.0]  # FIR filters have a_k = [1, 0, 0, ...]

    filters = {
        'Low-pass': (b_lp, a_lp),
        'High-pass': (b_hp, a_hp),
        'Band-pass': (b_bp, a_bp),
        'FIR': (b_fir, a_fir)
    }

    print(f"\nDigital Filter Design:")
    print(f"  Low-pass: Order {order_lp}, cutoff {cutoff_lp} Hz")
    print(f"  High-pass: Order {order_hp}, cutoff {cutoff_hp} Hz")
    print(f"  Band-pass: Order {order_bp}, band {low_band}-{high_band} Hz")
    print(f"  FIR: {numtaps} taps, cutoff {cutoff_fir} Hz")

    return filters

def apply_filters(signal_data, filters, fs):
    """Apply filters to signal and analyze results"""

    filtered_signals = {}

    for name, (b, a) in filters.items():
        # Apply filter
        filtered = signal.filtfilt(b, a, signal_data)
        filtered_signals[name] = filtered

        # Compute filter response
        w, h = signal.freqz(b, a, fs=fs)

        # Store frequency response
        filters[name] = (b, a, w, h)

    return filtered_signals

# Design filters and apply to test signal
filters = design_filters(fs)
filtered_sigs = apply_filters(test_sig, filters, fs)
\end{pycode}

\section{Window Functions and Spectral Estimation}
\label{sec:windows-spectral}

\begin{pycode}
# Window functions and spectral estimation
def compare_window_functions(N=256):
    """Compare different window functions"""

    window_funcs = {
        'Rectangular': np.ones(N),
        'Hamming': windows.hamming(N),
        'Hann': windows.hann(N),
        'Blackman': windows.blackman(N),
        'Kaiser': windows.kaiser(N, beta=8.6)
    }

    # Compute frequency responses
    window_spectra = {}
    freqs_win = np.linspace(-0.5, 0.5, 1024)

    for name, window in window_funcs.items():
        # Zero-pad for better frequency resolution
        padded = np.zeros(1024)
        padded[:N] = window

        # Compute FFT and shift
        spectrum = fftshift(fft(padded))
        window_spectra[name] = spectrum

        # Calculate metrics
        main_lobe_width = np.sum(np.abs(spectrum) > 0.5 * np.max(np.abs(spectrum))) / 1024
        side_lobe_level = 20 * np.log10(np.max(np.abs(spectrum[np.abs(spectrum) < 0.9*np.max(np.abs(spectrum))])) / np.max(np.abs(spectrum)))

        print(f"  {name}: Main lobe width ~= {main_lobe_width:.3f}, Side lobe level ~= {side_lobe_level:.1f} dB")

    print(f"\nWindow Function Analysis (N = {N}):")

    return window_funcs, window_spectra, freqs_win

def welch_spectral_estimation(signal_data, fs, nperseg=256, overlap_factor=0.5):
    """Welch's method for spectral estimation"""

    # Different window types
    window_types = ['hamming', 'hann', 'blackman']

    spectral_estimates = {}

    for window in window_types:
        freqs, psd = signal.welch(signal_data, fs, window=window,
                                nperseg=nperseg,
                                noverlap=int(nperseg * overlap_factor),
                                return_onesided=True)
        spectral_estimates[window] = (freqs, psd)

    print(f"\nWelch Spectral Estimation:")
    print(f"  Segment length: {nperseg}")
    print(f"  Overlap: {overlap_factor*100:.0f}%")
    print(f"  Windows tested: {window_types}")

    return spectral_estimates

# Analyze window functions and spectral estimation
window_funcs, win_spectra, freqs_win = compare_window_functions()
welch_estimates = welch_spectral_estimation(test_sig, fs)
\end{pycode}

\section{Time-Frequency Analysis}
\label{sec:time-frequency}

\begin{pycode}
# Time-frequency analysis using spectrograms and wavelets
def compute_spectrogram(signal_data, fs, nperseg=256, noverlap=None):
    """Compute and analyze spectrogram"""

    if noverlap is None:
        noverlap = nperseg // 2

    f, t_spec, Sxx = signal.spectrogram(signal_data, fs,
                                       nperseg=nperseg,
                                       noverlap=noverlap)

    print(f"\nSpectrogram Analysis:")
    print(f"  Time resolution: {t_spec[1] - t_spec[0]:.4f} s")
    print(f"  Frequency resolution: {f[1] - f[0]:.2f} Hz")
    print(f"  Time bins: {len(t_spec)}")
    print(f"  Frequency bins: {len(f)}")

    return f, t_spec, Sxx

def continuous_wavelet_transform(signal_data, fs, wavelet='morl'):
    """Simplified continuous wavelet transform demonstration"""

    # Define scales for CWT
    scales = np.arange(1, 128)

    # Create a simple Morlet-like wavelet manually
    def morlet_wavelet(M, w=5.0, s=1.0, complete=True):
        """
        Simple Morlet wavelet implementation
        """
        x = np.arange(0, M) - (M - 1.0) / 2
        x = x / s
        output = np.exp(1j * w * x) * np.exp(-0.5 * x**2) * np.pi**(-0.25)
        if not complete:
            output -= np.exp(-0.5 * w**2) * np.exp(-0.5 * x**2) * np.pi**(-0.25)
        return output

    # Compute CWT manually using convolution
    cwt_matrix = np.zeros((len(scales), len(signal_data)), dtype=complex)

    for i, scale in enumerate(scales):
        # Create wavelet at this scale
        wavelet_length = min(int(10 * scale), len(signal_data))
        if wavelet_length % 2 == 0:
            wavelet_length += 1

        psi = morlet_wavelet(wavelet_length, s=scale)

        # Normalize
        psi = psi / np.sqrt(scale)

        # Convolve signal with wavelet
        cwt_row = np.convolve(signal_data, np.conj(psi[::-1]), mode='same')
        cwt_matrix[i, :] = cwt_row

    # Convert scales to frequencies (approximate)
    frequencies = fs / (2 * scales)

    print(f"\nContinuous Wavelet Transform:")
    print(f"  Wavelet: {wavelet}")
    print(f"  Scales: {len(scales)}")
    print(f"  Frequency range: {frequencies[-1]:.1f} - {frequencies[0]:.1f} Hz")

    return frequencies, cwt_matrix

# Perform time-frequency analysis
f_spec, t_spec, Sxx = compute_spectrogram(test_sig, fs)
freq_cwt, cwt_coeffs = continuous_wavelet_transform(test_sig, fs)
\end{pycode}

\begin{pycode}
# Create comprehensive signal processing visualization
import os
os.makedirs('assets', exist_ok=True)

fig = plt.figure(figsize=(20, 15))
gs = fig.add_gridspec(4, 4, hspace=0.3, wspace=0.3)

# Time domain signal
ax1 = fig.add_subplot(gs[0, :2])
ax1.plot(t[:500], test_sig[:500], 'b-', linewidth=1.5, label='Test signal')
ax1.plot(t[:500], clean_sig[:500], 'r--', linewidth=1.5, alpha=0.7, label='Clean components')
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Amplitude')
ax1.set_title('Time Domain Signal')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Frequency domain (FFT)
ax2 = fig.add_subplot(gs[0, 2:])
ax2.loglog(pos_freqs[1:], psd_pos[1:], 'b-', linewidth=2)
ax2.set_xlabel('Frequency (Hz)')
ax2.set_ylabel('Power Spectral Density')
ax2.set_title('Power Spectrum (FFT)')
ax2.grid(True, alpha=0.3)

# Filter frequency responses
ax3 = fig.add_subplot(gs[1, :2])
for name, (b, a, w, h) in filters.items():
    if len(w) > 0:  # Check if frequency response was computed
        ax3.plot(w, 20*np.log10(np.abs(h)), linewidth=2, label=name)
ax3.set_xlabel('Frequency (Hz)')
ax3.set_ylabel('Magnitude (dB)')
ax3.set_title('Filter Frequency Responses')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Filtered signals (time domain)
ax4 = fig.add_subplot(gs[1, 2:])
time_slice = slice(0, 1000)  # Show first second
ax4.plot(t[time_slice], test_sig[time_slice], 'k-', alpha=0.5, label='Original')
colors = ['red', 'blue', 'green', 'orange']
for i, (name, filt_sig) in enumerate(filtered_sigs.items()):
    if i < len(colors):
        ax4.plot(t[time_slice], filt_sig[time_slice],
                color=colors[i], linewidth=1.5, label=name)
ax4.set_xlabel('Time (s)')
ax4.set_ylabel('Amplitude')
ax4.set_title('Filtered Signals')
ax4.legend()
ax4.grid(True, alpha=0.3)

# Window functions
ax5 = fig.add_subplot(gs[2, :2])
for name, window in window_funcs.items():
    ax5.plot(window, linewidth=2, label=name)
ax5.set_xlabel('Sample')
ax5.set_ylabel('Amplitude')
ax5.set_title('Window Functions')
ax5.legend()
ax5.grid(True, alpha=0.3)

# Window frequency responses
ax6 = fig.add_subplot(gs[2, 2:])
for name, spectrum in win_spectra.items():
    ax6.plot(freqs_win, 20*np.log10(np.abs(spectrum)/np.max(np.abs(spectrum))),
            linewidth=2, label=name)
ax6.set_xlabel('Normalized Frequency')
ax6.set_ylabel('Magnitude (dB)')
ax6.set_title('Window Frequency Responses')
ax6.set_ylim([-80, 5])
ax6.legend()
ax6.grid(True, alpha=0.3)

# Spectrogram
ax7 = fig.add_subplot(gs[3, :2])
pcm = ax7.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx), shading='gouraud', cmap='hot')
ax7.set_xlabel('Time (s)')
ax7.set_ylabel('Frequency (Hz)')
ax7.set_title('Spectrogram')
plt.colorbar(pcm, ax=ax7, label='Power (dB)')

# Wavelet transform
ax8 = fig.add_subplot(gs[3, 2:])
cwt_magnitude = np.abs(cwt_coeffs)
pcm2 = ax8.pcolormesh(t, freq_cwt, cwt_magnitude, shading='gouraud', cmap='viridis')
ax8.set_xlabel('Time (s)')
ax8.set_ylabel('Frequency (Hz)')
ax8.set_title('Continuous Wavelet Transform')
ax8.set_ylim([0, 400])  # Limit frequency range for visibility
plt.colorbar(pcm2, ax=ax8, label='|CWT|')

plt.suptitle('Signal Processing and Fourier Analysis', fontsize=16, fontweight='bold')
plt.savefig('assets/signal-processing-analysis.pdf', dpi=300, bbox_inches='tight')
plt.close()

print("Signal processing analysis saved to assets/signal-processing-analysis.pdf")
\end{pycode}

\begin{figure}[H]
    \centering
    \includegraphics[width=0.95\textwidth]{assets/signal-processing-analysis.pdf}
    \caption{Comprehensive signal processing and Fourier analysis including (top row) time domain signal with clean components and power spectrum via FFT, (second row) digital filter frequency responses and filtered signal comparisons, (third row) window function shapes and their frequency responses, and (bottom row) spectrogram and continuous wavelet transform providing time-frequency representations of the test signal.}
    \label{fig:signal-processing-analysis}
\end{figure}

\section{Conclusion}
\label{sec:conclusion}

This signal processing template demonstrates comprehensive spectral analysis techniques including:

\begin{itemize}
    \item Fourier transform theory and FFT implementation
    \item Digital filter design (IIR and FIR)
    \item Window functions and their spectral properties
    \item Time-frequency analysis via spectrograms and wavelets
    \item Professional visualization of signal characteristics
\end{itemize}

The computational methods provide a foundation for advanced signal processing applications in communications, audio processing, and scientific instrumentation.

\printbibliography

\end{document}