{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# SciPy Signal and Image Processing\n",
        "\n",
        "<!--\n",
        "Author: RSK World\n",
        "Website: https://rskworld.in\n",
        "Email: help@rskworld.in\n",
        "Phone: +91 93305 39277\n",
        "-->\n",
        "\n",
        "This notebook demonstrates signal and image processing with SciPy.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Author: RSK World - https://rskworld.in\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy import signal\n",
        "from scipy.ndimage import gaussian_filter, sobel\n",
        "from scipy.fft import fft, fftfreq\n",
        "\n",
        "print(\"SciPy Signal Processing Examples\")\n",
        "print(\"Author: RSK World - https://rskworld.in\")\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1. Signal Filtering\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2. Advanced: Wavelet Transform\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Author: RSK World - https://rskworld.in\n",
        "from scipy import signal\n",
        "\n",
        "# Generate signal with multiple time-frequency components\n",
        "np.random.seed(42)\n",
        "t = np.linspace(0, 1, 1000)\n",
        "signal_data = (np.sin(2 * np.pi * 10 * t) * (t < 0.5) + \n",
        "               np.sin(2 * np.pi * 25 * t) * (t >= 0.5) +\n",
        "               0.3 * np.random.randn(len(t)))\n",
        "\n",
        "# Continuous Wavelet Transform (using Morlet wavelet from scipy.signal)\n",
        "widths = np.arange(1, 31)\n",
        "cwt_matrix = signal.cwt(signal_data, signal.morlet, widths)\n",
        "\n",
        "# Visualize\n",
        "plt.figure(figsize=(14, 8))\n",
        "\n",
        "plt.subplot(3, 1, 1)\n",
        "plt.plot(t, signal_data, 'b-', linewidth=1.5)\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Amplitude', fontsize=12)\n",
        "plt.title('Original Signal', fontsize=12, fontweight='bold')\n",
        "plt.grid(True, alpha=0.3)\n",
        "\n",
        "plt.subplot(3, 1, 2)\n",
        "plt.imshow(np.abs(cwt_matrix), extent=[t[0], t[-1], widths[0], widths[-1]], \n",
        "           cmap='viridis', aspect='auto', origin='lower')\n",
        "plt.colorbar(label='Magnitude')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Scale', fontsize=12)\n",
        "plt.title('Continuous Wavelet Transform (CWT)', fontsize=12, fontweight='bold')\n",
        "\n",
        "# Short-Time Fourier Transform (STFT)\n",
        "f, t_stft, Zxx = signal.stft(signal_data, fs=1000, nperseg=100, noverlap=50)\n",
        "\n",
        "plt.subplot(3, 1, 3)\n",
        "plt.pcolormesh(t_stft, f, np.abs(Zxx), shading='gouraud', cmap='viridis')\n",
        "plt.colorbar(label='Magnitude')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Frequency (Hz)', fontsize=12)\n",
        "plt.title('Short-Time Fourier Transform (STFT)', fontsize=12, fontweight='bold')\n",
        "plt.ylim(0, 50)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3. Advanced: Image Processing - Morphological Operations\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Author: RSK World - https://rskworld.in\n",
        "from scipy.ndimage import binary_erosion, binary_dilation, binary_opening, binary_closing\n",
        "from scipy.ndimage import label, find_objects\n",
        "\n",
        "# Create test image with noise\n",
        "np.random.seed(42)\n",
        "image = np.zeros((100, 100))\n",
        "image[20:40, 20:40] = 1.0\n",
        "image[60:80, 60:80] = 1.0\n",
        "image_noisy = image + 0.3 * np.random.randn(100, 100)\n",
        "image_binary = (image_noisy > 0.5).astype(int)\n",
        "\n",
        "# Morphological operations\n",
        "eroded = binary_erosion(image_binary)\n",
        "dilated = binary_dilation(image_binary)\n",
        "opened = binary_opening(image_binary)\n",
        "closed = binary_closing(image_binary)\n",
        "\n",
        "# Connected components\n",
        "labeled, num_features = label(image_binary)\n",
        "\n",
        "# Visualize\n",
        "fig, axes = plt.subplots(2, 3, figsize=(15, 10))\n",
        "\n",
        "axes[0, 0].imshow(image_binary, cmap='gray')\n",
        "axes[0, 0].set_title('Binary Image', fontsize=12, fontweight='bold')\n",
        "axes[0, 0].axis('off')\n",
        "\n",
        "axes[0, 1].imshow(eroded, cmap='gray')\n",
        "axes[0, 1].set_title('Erosion', fontsize=12, fontweight='bold')\n",
        "axes[0, 1].axis('off')\n",
        "\n",
        "axes[0, 2].imshow(dilated, cmap='gray')\n",
        "axes[0, 2].set_title('Dilation', fontsize=12, fontweight='bold')\n",
        "axes[0, 2].axis('off')\n",
        "\n",
        "axes[1, 0].imshow(opened, cmap='gray')\n",
        "axes[1, 0].set_title('Opening (Erosion + Dilation)', fontsize=12, fontweight='bold')\n",
        "axes[1, 0].axis('off')\n",
        "\n",
        "axes[1, 1].imshow(closed, cmap='gray')\n",
        "axes[1, 1].set_title('Closing (Dilation + Erosion)', fontsize=12, fontweight='bold')\n",
        "axes[1, 1].axis('off')\n",
        "\n",
        "axes[1, 2].imshow(labeled, cmap='nipy_spectral')\n",
        "axes[1, 2].set_title(f'Connected Components (n={num_features})', fontsize=12, fontweight='bold')\n",
        "axes[1, 2].axis('off')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "print(f\"Number of connected components: {num_features}\")\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 4. Advanced: Spectral Analysis\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 5. Advanced: Image Denoising with Total Variation\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Author: RSK World - https://rskworld.in\n",
        "from scipy.ndimage import uniform_filter, gaussian_filter\n",
        "\n",
        "# Create test image\n",
        "np.random.seed(42)\n",
        "image_clean = np.zeros((100, 100))\n",
        "image_clean[20:40, 20:40] = 1.0\n",
        "image_clean[60:80, 60:80] = 1.0\n",
        "\n",
        "# Add heavy noise\n",
        "image_noisy = image_clean + 0.5 * np.random.randn(100, 100)\n",
        "image_noisy = np.clip(image_noisy, 0, 1)\n",
        "\n",
        "# Different denoising methods\n",
        "image_uniform = uniform_filter(image_noisy, size=5)\n",
        "image_gaussian = gaussian_filter(image_noisy, sigma=1.5)\n",
        "\n",
        "# Bilateral filter approximation (using combination)\n",
        "from scipy.ndimage import median_filter\n",
        "image_median = median_filter(image_noisy, size=3)\n",
        "\n",
        "# Visualize\n",
        "fig, axes = plt.subplots(2, 2, figsize=(12, 12))\n",
        "\n",
        "axes[0, 0].imshow(image_noisy, cmap='gray')\n",
        "axes[0, 0].set_title('Noisy Image', fontsize=12, fontweight='bold')\n",
        "axes[0, 0].axis('off')\n",
        "\n",
        "axes[0, 1].imshow(image_uniform, cmap='gray')\n",
        "axes[0, 1].set_title('Uniform Filter', fontsize=12, fontweight='bold')\n",
        "axes[0, 1].axis('off')\n",
        "\n",
        "axes[1, 0].imshow(image_gaussian, cmap='gray')\n",
        "axes[1, 0].set_title('Gaussian Filter', fontsize=12, fontweight='bold')\n",
        "axes[1, 0].axis('off')\n",
        "\n",
        "axes[1, 1].imshow(image_median, cmap='gray')\n",
        "axes[1, 1].set_title('Median Filter', fontsize=12, fontweight='bold')\n",
        "axes[1, 1].axis('off')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n",
        "# Calculate PSNR (Peak Signal-to-Noise Ratio)\n",
        "def psnr(img1, img2):\n",
        "    mse = np.mean((img1 - img2)**2)\n",
        "    if mse == 0:\n",
        "        return float('inf')\n",
        "    return 20 * np.log10(1.0 / np.sqrt(mse))\n",
        "\n",
        "print(f\"PSNR (Noisy): {psnr(image_clean, image_noisy):.2f} dB\")\n",
        "print(f\"PSNR (Uniform): {psnr(image_clean, image_uniform):.2f} dB\")\n",
        "print(f\"PSNR (Gaussian): {psnr(image_clean, image_gaussian):.2f} dB\")\n",
        "print(f\"PSNR (Median): {psnr(image_clean, image_median):.2f} dB\")\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 6. Advanced: Signal Decomposition (EMD-like)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Author: RSK World - https://rskworld.in\n",
        "from scipy.signal import hilbert\n",
        "\n",
        "# Generate composite signal\n",
        "np.random.seed(42)\n",
        "t = np.linspace(0, 2, 1000)\n",
        "fs = 1 / (t[1] - t[0])\n",
        "\n",
        "# Signal with multiple components\n",
        "signal_low = np.sin(2 * np.pi * 2 * t)  # Low frequency\n",
        "signal_high = 0.5 * np.sin(2 * np.pi * 10 * t)  # High frequency\n",
        "signal_trend = 0.1 * t  # Trend\n",
        "noise = 0.1 * np.random.randn(len(t))\n",
        "\n",
        "signal_composite = signal_low + signal_high + signal_trend + noise\n",
        "\n",
        "# Extract components using filtering\n",
        "b_low, a_low = signal.butter(4, 0.05, 'low')\n",
        "signal_low_extracted = signal.filtfilt(b_low, a_low, signal_composite)\n",
        "\n",
        "b_high, a_high = signal.butter(4, [0.08, 0.15], 'band')\n",
        "signal_band_extracted = signal.filtfilt(b_high, a_high, signal_composite)\n",
        "\n",
        "# Trend extraction (high-pass filter)\n",
        "b_trend, a_trend = signal.butter(4, 0.01, 'high')\n",
        "signal_trend_extracted = signal.filtfilt(b_trend, a_trend, signal_composite)\n",
        "\n",
        "# Hilbert transform for instantaneous frequency\n",
        "analytic_signal = hilbert(signal_low_extracted)\n",
        "instantaneous_phase = np.unwrap(np.angle(analytic_signal))\n",
        "instantaneous_frequency = (np.diff(instantaneous_phase) / (2.0 * np.pi) * fs)\n",
        "\n",
        "# Visualize\n",
        "plt.figure(figsize=(14, 10))\n",
        "\n",
        "plt.subplot(4, 1, 1)\n",
        "plt.plot(t, signal_composite, 'b-', linewidth=1.5)\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Amplitude', fontsize=12)\n",
        "plt.title('Composite Signal', fontsize=12, fontweight='bold')\n",
        "plt.grid(True, alpha=0.3)\n",
        "\n",
        "plt.subplot(4, 1, 2)\n",
        "plt.plot(t, signal_low_extracted, 'r-', linewidth=1.5, label='Low frequency')\n",
        "plt.plot(t, signal_low, 'r--', linewidth=1, alpha=0.5, label='True low')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Amplitude', fontsize=12)\n",
        "plt.title('Low Frequency Component', fontsize=12, fontweight='bold')\n",
        "plt.legend()\n",
        "plt.grid(True, alpha=0.3)\n",
        "\n",
        "plt.subplot(4, 1, 3)\n",
        "plt.plot(t, signal_band_extracted, 'g-', linewidth=1.5, label='Band frequency')\n",
        "plt.plot(t, signal_high, 'g--', linewidth=1, alpha=0.5, label='True high')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Amplitude', fontsize=12)\n",
        "plt.title('High Frequency Component', fontsize=12, fontweight='bold')\n",
        "plt.legend()\n",
        "plt.grid(True, alpha=0.3)\n",
        "\n",
        "plt.subplot(4, 1, 4)\n",
        "plt.plot(t[1:], instantaneous_frequency, 'm-', linewidth=1.5)\n",
        "plt.axhline(2, color='k', linestyle='--', linewidth=1, label='True frequency')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Frequency (Hz)', fontsize=12)\n",
        "plt.title('Instantaneous Frequency', fontsize=12, fontweight='bold')\n",
        "plt.legend()\n",
        "plt.grid(True, alpha=0.3)\n",
        "plt.ylim(0, 5)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Advanced spectral analysis\n",
        "from scipy.signal import welch, spectrogram\n",
        "\n",
        "# Generate signal with multiple frequencies\n",
        "np.random.seed(42)\n",
        "t = np.linspace(0, 2, 2000)\n",
        "fs = 1 / (t[1] - t[0])\n",
        "signal_data = (np.sin(2 * np.pi * 10 * t) + \n",
        "               0.5 * np.sin(2 * np.pi * 25 * t) +\n",
        "               0.3 * np.sin(2 * np.pi * 50 * t) +\n",
        "               0.2 * np.random.randn(len(t)))\n",
        "\n",
        "# Power Spectral Density using Welch's method\n",
        "f_psd, psd = welch(signal_data, fs, nperseg=256)\n",
        "\n",
        "# Spectrogram\n",
        "f_spec, t_spec, Sxx = spectrogram(signal_data, fs, nperseg=128, noverlap=64)\n",
        "\n",
        "# Visualize\n",
        "plt.figure(figsize=(14, 8))\n",
        "\n",
        "plt.subplot(3, 1, 1)\n",
        "plt.plot(t, signal_data, 'b-', linewidth=1)\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Amplitude', fontsize=12)\n",
        "plt.title('Time Domain Signal', fontsize=12, fontweight='bold')\n",
        "plt.grid(True, alpha=0.3)\n",
        "\n",
        "plt.subplot(3, 1, 2)\n",
        "plt.semilogy(f_psd, psd, 'r-', linewidth=2)\n",
        "plt.xlabel('Frequency (Hz)', fontsize=12)\n",
        "plt.ylabel('Power Spectral Density', fontsize=12)\n",
        "plt.title('Power Spectral Density (Welch Method)', fontsize=12, fontweight='bold')\n",
        "plt.grid(True, alpha=0.3)\n",
        "plt.xlim(0, 60)\n",
        "\n",
        "plt.subplot(3, 1, 3)\n",
        "plt.pcolormesh(t_spec, f_spec, 10 * np.log10(Sxx), shading='gouraud', cmap='viridis')\n",
        "plt.colorbar(label='Power (dB)')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Frequency (Hz)', fontsize=12)\n",
        "plt.title('Spectrogram', fontsize=12, fontweight='bold')\n",
        "plt.ylim(0, 60)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Generate noisy signal\n",
        "np.random.seed(42)\n",
        "t = np.linspace(0, 1, 1000)\n",
        "signal_clean = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 10 * t)\n",
        "noise = 0.5 * np.random.randn(len(t))\n",
        "signal_noisy = signal_clean + noise\n",
        "\n",
        "# Low-pass filter\n",
        "b, a = signal.butter(4, 0.1, 'low')\n",
        "signal_filtered = signal.filtfilt(b, a, signal_noisy)\n",
        "\n",
        "# Visualize\n",
        "plt.figure(figsize=(14, 6))\n",
        "\n",
        "plt.subplot(2, 1, 1)\n",
        "plt.plot(t, signal_clean, 'b-', linewidth=2, label='Clean signal')\n",
        "plt.plot(t, signal_noisy, 'r-', linewidth=1, alpha=0.5, label='Noisy signal')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Amplitude', fontsize=12)\n",
        "plt.title('Original Signals', fontsize=12, fontweight='bold')\n",
        "plt.legend()\n",
        "plt.grid(True, alpha=0.3)\n",
        "\n",
        "plt.subplot(2, 1, 2)\n",
        "plt.plot(t, signal_clean, 'b--', linewidth=2, alpha=0.7, label='Clean signal')\n",
        "plt.plot(t, signal_filtered, 'g-', linewidth=2, label='Filtered signal')\n",
        "plt.xlabel('Time (s)', fontsize=12)\n",
        "plt.ylabel('Amplitude', fontsize=12)\n",
        "plt.title('Butterworth Low-Pass Filter', fontsize=12, fontweight='bold')\n",
        "plt.legend()\n",
        "plt.grid(True, alpha=0.3)\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()\n",
        "\n"
      ]
    }
  ],
  "metadata": {
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}