"""
Signal and Image Processing with SciPy
Author: RSK World
Website: https://rskworld.in
Email: help@rskworld.in
Phone: +91 93305 39277
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.ndimage import gaussian_filter, sobel, median_filter
from scipy.fft import fft, fftfreq, ifft

# Example 1: Signal Filtering
def example_signal_filtering():
    """
    Signal filtering examples
    """
    print("=" * 60)
    print("Example 1: Signal Filtering")
    print("=" * 60)
    
    # Generate noisy signal
    np.random.seed(42)
    t = np.linspace(0, 1, 1000)
    signal_clean = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 10 * t)
    noise = 0.5 * np.random.randn(len(t))
    signal_noisy = signal_clean + noise
    
    # Low-pass filter
    b, a = signal.butter(4, 0.1, 'low')
    signal_filtered = signal.filtfilt(b, a, signal_noisy)
    
    # Moving average filter
    window_size = 20
    signal_ma = np.convolve(signal_noisy, np.ones(window_size)/window_size, mode='same')
    
    print(f"Signal length: {len(t)}")
    print(f"Sampling frequency: {1/(t[1]-t[0]):.2f} Hz")
    print(f"Filter types: Butterworth low-pass, Moving average")
    
    # Visualize
    plt.figure(figsize=(14, 8))
    
    plt.subplot(3, 1, 1)
    plt.plot(t, signal_clean, 'b-', linewidth=2, label='Clean signal')
    plt.plot(t, signal_noisy, 'r-', linewidth=1, alpha=0.5, label='Noisy signal')
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Original Signals', fontsize=12, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.subplot(3, 1, 2)
    plt.plot(t, signal_clean, 'b--', linewidth=2, alpha=0.7, label='Clean signal')
    plt.plot(t, signal_filtered, 'g-', linewidth=2, label='Butterworth filter')
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Butterworth Low-Pass Filter', fontsize=12, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.subplot(3, 1, 3)
    plt.plot(t, signal_clean, 'b--', linewidth=2, alpha=0.7, label='Clean signal')
    plt.plot(t, signal_ma, 'm-', linewidth=2, label='Moving average filter')
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Moving Average Filter', fontsize=12, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('signal_filtering.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'signal_filtering.png'")
    plt.close()


# Example 2: Frequency Domain Analysis
def example_frequency_analysis():
    """
    Frequency domain analysis using FFT
    """
    print("\n" + "=" * 60)
    print("Example 2: Frequency Domain Analysis")
    print("=" * 60)
    
    # Generate signal with multiple frequencies
    np.random.seed(42)
    t = np.linspace(0, 1, 1000)
    dt = t[1] - t[0]
    fs = 1 / dt
    
    # Signal: sum of sinusoids at 5 Hz, 10 Hz, and 20 Hz
    signal_time = (np.sin(2 * np.pi * 5 * t) + 
                   0.5 * np.sin(2 * np.pi * 10 * t) + 
                   0.3 * np.sin(2 * np.pi * 20 * t) +
                   0.2 * np.random.randn(len(t)))
    
    # FFT
    fft_vals = fft(signal_time)
    fft_freq = fftfreq(len(signal_time), dt)
    
    # Power spectral density
    psd = np.abs(fft_vals)**2
    
    # Only positive frequencies
    positive_freq_idx = fft_freq > 0
    fft_freq_pos = fft_freq[positive_freq_idx]
    psd_pos = psd[positive_freq_idx]
    
    print(f"Sampling frequency: {fs:.2f} Hz")
    print(f"Signal frequencies: 5 Hz, 10 Hz, 20 Hz")
    print(f"Peak frequencies detected:")
    peak_indices = signal.find_peaks(psd_pos, height=np.max(psd_pos)*0.1)[0]
    for idx in peak_indices[:5]:
        print(f"  {fft_freq_pos[idx]:.2f} Hz")
    
    # Visualize
    plt.figure(figsize=(14, 6))
    
    plt.subplot(1, 2, 1)
    plt.plot(t, signal_time, 'b-', linewidth=1.5)
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Time Domain Signal', fontsize=12, fontweight='bold')
    plt.grid(True, alpha=0.3)
    
    plt.subplot(1, 2, 2)
    plt.plot(fft_freq_pos, psd_pos, 'r-', linewidth=2)
    plt.xlabel('Frequency (Hz)', fontsize=12)
    plt.ylabel('Power Spectral Density', fontsize=12)
    plt.title('Frequency Domain (FFT)', fontsize=12, fontweight='bold')
    plt.xlim(0, 50)
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('signal_frequency.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'signal_frequency.png'")
    plt.close()


# Example 3: Image Filtering
def example_image_filtering():
    """
    Image filtering examples
    """
    print("\n" + "=" * 60)
    print("Example 3: Image Filtering")
    print("=" * 60)
    
    # Create a test image with noise
    np.random.seed(42)
    image_clean = np.zeros((100, 100))
    image_clean[20:40, 20:40] = 1.0
    image_clean[60:80, 60:80] = 1.0
    
    # Add noise
    image_noisy = image_clean + 0.3 * np.random.randn(100, 100)
    image_noisy = np.clip(image_noisy, 0, 1)
    
    # Apply filters
    image_gaussian = gaussian_filter(image_noisy, sigma=1.0)
    image_median = median_filter(image_noisy, size=3)
    
    print(f"Image size: {image_clean.shape}")
    print(f"Filters applied: Gaussian, Median")
    
    # Visualize
    fig, axes = plt.subplots(2, 2, figsize=(12, 12))
    
    axes[0, 0].imshow(image_clean, cmap='gray')
    axes[0, 0].set_title('Clean Image', fontsize=12, fontweight='bold')
    axes[0, 0].axis('off')
    
    axes[0, 1].imshow(image_noisy, cmap='gray')
    axes[0, 1].set_title('Noisy Image', fontsize=12, fontweight='bold')
    axes[0, 1].axis('off')
    
    axes[1, 0].imshow(image_gaussian, cmap='gray')
    axes[1, 0].set_title('Gaussian Filter', fontsize=12, fontweight='bold')
    axes[1, 0].axis('off')
    
    axes[1, 1].imshow(image_median, cmap='gray')
    axes[1, 1].set_title('Median Filter', fontsize=12, fontweight='bold')
    axes[1, 1].axis('off')
    
    plt.tight_layout()
    plt.savefig('signal_image_filtering.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'signal_image_filtering.png'")
    plt.close()


# Example 4: Edge Detection
def example_edge_detection():
    """
    Edge detection in images
    """
    print("\n" + "=" * 60)
    print("Example 4: Edge Detection")
    print("=" * 60)
    
    # Create test image with edges
    image = np.zeros((100, 100))
    image[30:70, 30:70] = 1.0
    image[10:20, 10:90] = 0.5
    
    # Add noise
    np.random.seed(42)
    image_noisy = image + 0.1 * np.random.randn(100, 100)
    image_noisy = np.clip(image_noisy, 0, 1)
    
    # Edge detection using Sobel operator
    edges_x = sobel(image_noisy, axis=1)
    edges_y = sobel(image_noisy, axis=0)
    edges_magnitude = np.sqrt(edges_x**2 + edges_y**2)
    
    print(f"Edge detection method: Sobel operator")
    print(f"Image size: {image_noisy.shape}")
    
    # Visualize
    fig, axes = plt.subplots(2, 2, figsize=(12, 12))
    
    axes[0, 0].imshow(image_noisy, cmap='gray')
    axes[0, 0].set_title('Original Image', fontsize=12, fontweight='bold')
    axes[0, 0].axis('off')
    
    axes[0, 1].imshow(edges_x, cmap='gray')
    axes[0, 1].set_title('Edges (X-direction)', fontsize=12, fontweight='bold')
    axes[0, 1].axis('off')
    
    axes[1, 0].imshow(edges_y, cmap='gray')
    axes[1, 0].set_title('Edges (Y-direction)', fontsize=12, fontweight='bold')
    axes[1, 0].axis('off')
    
    axes[1, 1].imshow(edges_magnitude, cmap='hot')
    axes[1, 1].set_title('Edge Magnitude', fontsize=12, fontweight='bold')
    axes[1, 1].axis('off')
    
    plt.tight_layout()
    plt.savefig('signal_edge_detection.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'signal_edge_detection.png'")
    plt.close()


# Example 5: Signal Convolution
def example_convolution():
    """
    Signal convolution examples
    """
    print("\n" + "=" * 60)
    print("Example 5: Signal Convolution")
    print("=" * 60)
    
    # Create signals
    t = np.linspace(0, 5, 500)
    
    # Input signal: rectangular pulse
    signal_input = np.zeros_like(t)
    signal_input[(t >= 1) & (t <= 2)] = 1.0
    
    # Impulse response: exponential decay
    h = np.exp(-3 * t)
    h = h / np.sum(h)  # Normalize
    
    # Convolution
    signal_output = np.convolve(signal_input, h, mode='same')
    
    print(f"Input signal: Rectangular pulse")
    print(f"Impulse response: Exponential decay")
    print(f"Output: Convolution of input and impulse response")
    
    # Visualize
    plt.figure(figsize=(14, 8))
    
    plt.subplot(3, 1, 1)
    plt.plot(t, signal_input, 'b-', linewidth=2)
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Input Signal', fontsize=12, fontweight='bold')
    plt.grid(True, alpha=0.3)
    plt.ylim(-0.1, 1.2)
    
    plt.subplot(3, 1, 2)
    plt.plot(t, h, 'g-', linewidth=2)
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Impulse Response', fontsize=12, fontweight='bold')
    plt.grid(True, alpha=0.3)
    
    plt.subplot(3, 1, 3)
    plt.plot(t, signal_output, 'r-', linewidth=2)
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Output Signal (Convolution)', fontsize=12, fontweight='bold')
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('signal_convolution.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'signal_convolution.png'")
    plt.close()


# Example 6: Signal Detection
def example_signal_detection():
    """
    Signal detection using correlation
    """
    print("\n" + "=" * 60)
    print("Example 6: Signal Detection")
    print("=" * 60)
    
    # Generate signal with embedded pattern
    np.random.seed(42)
    t = np.linspace(0, 10, 1000)
    
    # Pattern to detect: sinusoidal burst
    pattern = np.sin(2 * np.pi * 5 * np.linspace(0, 0.5, 50))
    
    # Create signal with pattern at specific location
    signal_full = 0.1 * np.random.randn(len(t))
    pattern_start = 300
    signal_full[pattern_start:pattern_start+len(pattern)] += pattern
    
    # Cross-correlation for detection
    correlation = np.correlate(signal_full, pattern, mode='same')
    correlation_normalized = correlation / np.max(np.abs(correlation))
    
    # Find peak
    peak_idx = np.argmax(np.abs(correlation_normalized))
    peak_time = t[peak_idx]
    
    print(f"Pattern length: {len(pattern)} samples")
    print(f"True pattern location: {t[pattern_start]:.2f} s")
    print(f"Detected location: {peak_time:.2f} s")
    print(f"Detection error: {abs(peak_time - t[pattern_start]):.4f} s")
    
    # Visualize
    plt.figure(figsize=(14, 8))
    
    plt.subplot(3, 1, 1)
    plt.plot(t, signal_full, 'b-', linewidth=1)
    plt.axvline(t[pattern_start], color='red', linestyle='--', linewidth=2, label='True location')
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Signal with Embedded Pattern', fontsize=12, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.subplot(3, 1, 2)
    pattern_t = np.linspace(0, 0.5, len(pattern))
    plt.plot(pattern_t, pattern, 'g-', linewidth=2)
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Amplitude', fontsize=12)
    plt.title('Pattern to Detect', fontsize=12, fontweight='bold')
    plt.grid(True, alpha=0.3)
    
    plt.subplot(3, 1, 3)
    plt.plot(t, correlation_normalized, 'r-', linewidth=2)
    plt.axvline(peak_time, color='red', linestyle='--', linewidth=2, label='Detected location')
    plt.xlabel('Time (s)', fontsize=12)
    plt.ylabel('Correlation', fontsize=12)
    plt.title('Cross-Correlation (Detection)', fontsize=12, fontweight='bold')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('signal_detection.png', dpi=300, bbox_inches='tight')
    print("\nPlot saved as 'signal_detection.png'")
    plt.close()


def main():
    """
    Main function to run all signal processing examples
    """
    print("\n" + "=" * 60)
    print("SciPy Signal and Image Processing Examples")
    print("Author: RSK World - https://rskworld.in")
    print("=" * 60 + "\n")
    
    example_signal_filtering()
    example_frequency_analysis()
    example_image_filtering()
    example_edge_detection()
    example_convolution()
    example_signal_detection()
    
    print("\n" + "=" * 60)
    print("All signal processing examples completed!")
    print("=" * 60)


if __name__ == "__main__":
    main()