ubuntu2404
=>PYTHONTEX#py#default#default#0#code#####42#
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)
=>PYTHONTEX#py#default#default#1#code#####134#
# 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)
=>PYTHONTEX#py#default#default#2#code#####203#
# 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)
=>PYTHONTEX#py#default#default#3#code#####269#
# 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)
=>PYTHONTEX#py#default#default#4#code#####340#
# 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")
=>PYTHONTEX:SETTINGS#
version=0.18
outputdir=pythontex-files-main
workingdir=.
workingdirset=false
gobble=none
rerun=default
hashdependencies=default
makestderr=false
stderrfilename=full
keeptemps=none
pyfuture=default
pyconfuture=none
pygments=true
pygglobal=:GLOBAL||
fvextfile=-1
pyconbanner=none
pyconfilename=stdin
depythontex=false
pygfamily=py|python3|
pygfamily=pycon|pycon|
pygfamily=sympy|python3|
pygfamily=sympycon|pycon|
pygfamily=pylab|python3|
pygfamily=pylabcon|pycon|