{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ed23fea8",
   "metadata": {},
   "source": [
    "# Signals and Systems (ECE-Y425)\n",
    "\n",
    "## Electrical and Computer Engineering Department\n",
    "\n",
    "### University of Patras, Greece\n",
    "\n",
    "**Instructor:** Konstantinos Chatzilygeroudis\n",
    "\n",
    "**Email:** [costashatz@upatras.gr](mailto:costashatz@upatras.gr)\n",
    "\n",
    "#### Hands-On Session on Fourier Analysis on Realistic Signals\n",
    "\n",
    "**Objective:** We are given a signal that is changing over time, and we want to identify its frequency components."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fef83dfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (18, 4)\n",
    "plt.rcParams[\"axes.grid\"] = True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5ec9829",
   "metadata": {},
   "source": [
    "## 1. Build a realistic non-stationary signal\n",
    "\n",
    "We construct a 3-second signal:\n",
    "\n",
    "- Segment 1: 100 Hz tone\n",
    "- Segment 2: 100 Hz + 180 Hz\n",
    "- Segment 3: linear chirp from 100 Hz to 300 Hz\n",
    "- Additive white noise everywhere\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db5084a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sampling parameters\n",
    "fs = 2000            # Hz\n",
    "T_total = 3.0        # seconds\n",
    "N = int(fs * T_total)\n",
    "t = np.arange(N) / fs\n",
    "\n",
    "# Segment masks\n",
    "seg1 = (t < 1.0)\n",
    "seg2 = (t >= 1.0) & (t < 2.0)\n",
    "seg3 = (t >= 2.0)\n",
    "\n",
    "# Create empty signal\n",
    "x = np.zeros_like(t)\n",
    "\n",
    "# Segment 1: one tone\n",
    "x[seg1] = 1.0 * np.sin(2*np.pi*100*t[seg1])\n",
    "\n",
    "# Segment 2: two tones\n",
    "x[seg2] = (\n",
    "    0.9 * np.sin(2*np.pi*100*t[seg2]) +\n",
    "    0.7 * np.sin(2*np.pi*180*t[seg2])\n",
    ")\n",
    "\n",
    "# Segment 3: chirp 100 -> 300 Hz over 1 second\n",
    "t3 = t[seg3] - 2.0\n",
    "f0 = 100\n",
    "f1 = 300\n",
    "T3 = 1.0\n",
    "k = (f1 - f0) / T3     # chirp rate in Hz/s\n",
    "phase3 = 2*np.pi*(f0*t3 + 0.5*k*t3**2)\n",
    "x[seg3] = 0.9 * np.sin(phase3)\n",
    "\n",
    "# Add noise\n",
    "rng = np.random.default_rng(7)\n",
    "noise_std = 0.35\n",
    "x_noisy = x + noise_std * rng.standard_normal(N)\n",
    "\n",
    "print(f\"Signal length: {T_total:.1f} s\")\n",
    "print(f\"Sampling rate: {fs} Hz\")\n",
    "print(f\"Number of samples: {N}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78bb8a4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(t, x_noisy, label=\"Noisy signal\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.title(\"Time-domain signal\")\n",
    "# plt.xlim([0., 0.3])\n",
    "# plt.xlim([1., 1.3])\n",
    "# plt.xlim([2., 2.3])\n",
    "# plt.xlim([2.7, 3.])\n",
    "# plt.xlim([2., 3.])\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caa6e7a8",
   "metadata": {},
   "source": [
    "### Inspect the time waveform.\n",
    "\n",
    "- Can you clearly identify all three behaviors from the time plot alone?\n",
    "- Where does the signal seem to change?\n",
    "- Why is this difficult?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0af39912",
   "metadata": {},
   "source": [
    "## 2. Useful helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfc74903",
   "metadata": {},
   "outputs": [],
   "source": [
    "def one_sided_fft(x, fs):\n",
    "    \"\"\"Return one-sided frequency axis and magnitude spectrum.\"\"\"\n",
    "    N = len(x)\n",
    "    X = np.fft.rfft(x)\n",
    "    f = np.fft.rfftfreq(N, d=1/fs)\n",
    "    mag = np.abs(X)\n",
    "    return f, mag\n",
    "\n",
    "def mag_to_db(mag, eps=1e-12):\n",
    "    return 20*np.log10(np.maximum(mag, eps))\n",
    "\n",
    "def coherent_gain(w):\n",
    "    return np.mean(w)\n",
    "\n",
    "def apply_window(x, window_name=\"rect\"):\n",
    "    N = len(x)\n",
    "    if window_name == \"rect\":\n",
    "        w = np.ones(N)\n",
    "    elif window_name == \"hann\":\n",
    "        w = np.hanning(N)\n",
    "    elif window_name == \"hamming\":\n",
    "        w = np.hamming(N)\n",
    "    elif window_name == \"blackman\":\n",
    "        w = np.blackman(N)\n",
    "    else:\n",
    "        raise ValueError(\"Unknown window.\")\n",
    "    return x * w, w\n",
    "\n",
    "def plot_spectrum(f, mag, title=\"\", xlim=None, db=True, normalize_db=False):\n",
    "    plt.figure()\n",
    "    if db:\n",
    "        if normalize_db:\n",
    "            mg = mag_to_db(mag / np.max(np.abs(mag)))\n",
    "        else:\n",
    "            mg = mag_to_db(mag)\n",
    "        plt.plot(f, mg)\n",
    "        plt.ylabel(\"Magnitude (dB)\")\n",
    "    else:\n",
    "        plt.plot(f, mag)\n",
    "        plt.ylabel(\"Magnitude\")\n",
    "    plt.xlabel(\"Frequency (Hz)\")\n",
    "    plt.title(title)\n",
    "    if xlim is not None:\n",
    "        plt.xlim(xlim)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89fd6ba5",
   "metadata": {},
   "source": [
    "## 3. Full-record FFT: the naive analysis\n",
    "\n",
    "A first instinct is to compute the FFT of the **entire 3-second signal**.\n",
    "\n",
    "This is useful — but because the signal changes over time, the FFT mixes all behaviors together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fec2502b",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_full, mag_full = one_sided_fft(x_noisy, fs)\n",
    "plot_spectrum(f_full, mag_full, title=\"FFT of the full 3-second signal\", xlim=(0, 400))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a5e4c2b",
   "metadata": {},
   "source": [
    "### Inspect the spectrum\n",
    "\n",
    "- Which frequencies can you identify confidently?\n",
    "- Can you tell **when** each component exists?\n",
    "- Can you distinguish the chirp from a set of stationary tones?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "108c2cf5",
   "metadata": {},
   "source": [
    "## 4. Segment-based analysis (manual short-time Fourier analysis)\n",
    "\n",
    "To reveal time variation, we split the signal into short chunks and compute an FFT for each chunk.\n",
    "\n",
    "This is the core idea behind a **spectrogram**, but here we do it manually first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97105ce6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def segment_signal(x, fs, segment_duration, hop_duration=None):\n",
    "    \"\"\"Return segments and their start times.\"\"\"\n",
    "    if hop_duration is None:\n",
    "        hop_duration = segment_duration\n",
    "\n",
    "    L = int(segment_duration * fs)\n",
    "    H = int(hop_duration * fs)\n",
    "\n",
    "    segments = []\n",
    "    starts = []\n",
    "    for start in range(0, len(x) - L + 1, H):\n",
    "        segments.append(x[start:start+L])\n",
    "        starts.append(start / fs)\n",
    "    return np.array(segments), np.array(starts)\n",
    "\n",
    "segment_duration = 0.25   # 250 ms\n",
    "hop_duration = 0.125      # 125 ms\n",
    "\n",
    "segments, starts = segment_signal(x_noisy, fs, segment_duration, hop_duration)\n",
    "\n",
    "print(\"Number of segments:\", len(segments))\n",
    "print(\"Segment length (samples):\", segments.shape[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb7be227",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize the segments in time\n",
    "plt.figure(figsize=(18, 6))\n",
    "for i in range(len(segments)):\n",
    "    plt.plot(np.arange(len(segments[i]))/fs + starts[i], segments[i], label=f\"Seg {i}\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "plt.title(\"Example segments\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e231413",
   "metadata": {},
   "source": [
    "## 5. Analyze one segment at a time\n",
    "\n",
    "We start with a rectangular window, then later compare other windows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ead2764c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose a few representative segments\n",
    "chosen_indices = [1, 8, 18]   # roughly from the three regimes\n",
    "\n",
    "for idx in chosen_indices:\n",
    "    xseg = segments[idx]\n",
    "    f, mag = one_sided_fft(xseg, fs)\n",
    "    plot_spectrum(\n",
    "        f, mag,\n",
    "        title=f\"Segment starting at t = {starts[idx]:.3f} s (rectangular window)\",\n",
    "        xlim=(0, 400)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68fe7971",
   "metadata": {},
   "source": [
    "### Estimating the dominant frequency for each segment\n",
    "For each of the three selected segments:\n",
    "\n",
    "- Estimate the dominant frequencies\n",
    "- Explain how the spectra differ\n",
    "- Which segment corresponds to the chirp?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b940af7b",
   "metadata": {},
   "source": [
    "## 6. Track the dominant frequency over time\n",
    "\n",
    "A simple way to summarize each segment is to find the frequency of the largest FFT peak.\n",
    "This is crude, but very intuitive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8199773c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dominant_frequency(xseg, fs, window_name=\"rect\", fmax=None):\n",
    "    xw, w = apply_window(xseg, window_name)\n",
    "    f, mag = one_sided_fft(xw, fs)\n",
    "    cg = coherent_gain(w)\n",
    "    mag = mag / max(cg, 1e-12)  # crude amplitude correction\n",
    "\n",
    "    if fmax is not None:\n",
    "        valid = f <= fmax\n",
    "        f = f[valid]\n",
    "        mag = mag[valid]\n",
    "\n",
    "    idx = np.argmax(mag)\n",
    "    return f[idx], mag[idx]\n",
    "\n",
    "dom_freqs = []\n",
    "for xseg in segments:\n",
    "    fd, _ = dominant_frequency(xseg, fs, window_name=\"hann\", fmax=400)\n",
    "    dom_freqs.append(fd)\n",
    "\n",
    "dom_freqs = np.array(dom_freqs)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(starts + segment_duration/2, dom_freqs, marker=\"o\")\n",
    "plt.xlabel(\"Segment center time (s)\")\n",
    "plt.ylabel(\"Dominant frequency (Hz)\")\n",
    "plt.title(\"Dominant frequency vs time\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34ced014",
   "metadata": {},
   "source": [
    "### Interpret the dominant-frequency plot.\n",
    "\n",
    "- Where do you see transitions?\n",
    "- Why does the “two-tone” region still collapse to one dominant frequency?\n",
    "- What limitations does this method have?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3e9b542",
   "metadata": {},
   "source": [
    "## 7. Window comparison on the same segment\n",
    "\n",
    "Now let's compare different windows on one segment.\n",
    "\n",
    "We pick a segment where leakage matters — for example in the two-tone region or chirp region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93adf847",
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = 8   # try changing this\n",
    "xseg = segments[idx]\n",
    "\n",
    "for win_name in [\"rect\", \"hann\", \"hamming\", \"blackman\"]:\n",
    "    xw, w = apply_window(xseg, win_name)\n",
    "    f, mag = one_sided_fft(xw, fs)\n",
    "    mag = mag / coherent_gain(w)\n",
    "    plot_spectrum(\n",
    "        f, mag,\n",
    "        title=f\"Window = {win_name}, segment starting at {starts[idx]:.3f} s\",\n",
    "        xlim=(0, 400)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2196aa1a",
   "metadata": {},
   "source": [
    "### Compare the windows:\n",
    "\n",
    "- Which one gives the cleanest sidelobe behavior?\n",
    "- Which one makes peaks look wider?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbc4efe0",
   "metadata": {},
   "source": [
    "## 8. Time–frequency tradeoff\n",
    "\n",
    "Now, we try different segment lengths. Short windows give better time localization; long windows give better frequency resolution.\n",
    "\n",
    "Let's compare:\n",
    "\n",
    "- 50 ms\n",
    "- 100 ms\n",
    "- 500 ms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd0a0e03",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dominant_frequency_curve(x, fs, segment_duration, hop_duration, window_name=\"hann\", fmax=400):\n",
    "    segs, sts = segment_signal(x, fs, segment_duration, hop_duration)\n",
    "    dom = []\n",
    "    for s in segs:\n",
    "        fd, _ = dominant_frequency(s, fs, window_name=window_name, fmax=fmax)\n",
    "        dom.append(fd)\n",
    "    return sts + segment_duration/2, np.array(dom)\n",
    "\n",
    "configs = [\n",
    "    (0.05, 0.025),\n",
    "    (0.10, 0.05),\n",
    "    (0.50, 0.25),\n",
    "]\n",
    "\n",
    "plt.figure()\n",
    "for seg_dur, hop_dur in configs:\n",
    "    tc, df = dominant_frequency_curve(x_noisy, fs, seg_dur, hop_dur, window_name=\"hann\", fmax=400)\n",
    "    plt.plot(tc, df, marker=\"o\", label=f\"{int(seg_dur*1000)} ms\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"Dominant frequency (Hz)\")\n",
    "plt.title(\"Effect of segment duration\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b32c732",
   "metadata": {},
   "source": [
    "### Explain the tradeoff:\n",
    "\n",
    "- Which segment duration tracks the transitions most sharply?\n",
    "- Which segment duration gives the most stable frequency estimate?\n",
    "- Why can’t we get both perfectly at the same time?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36ca74c3",
   "metadata": {},
   "source": [
    "## 9. Build a manual spectral map\n",
    "\n",
    "Instead of keeping only the largest peak, compute the full FFT magnitude for each segment and stack the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "346524ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def spectral_map(x, fs, segment_duration, hop_duration, window_name=\"hann\", fmax=400):\n",
    "    segs, sts = segment_signal(x, fs, segment_duration, hop_duration)\n",
    "    rows = []\n",
    "    freq_axis = None\n",
    "\n",
    "    for s in segs:\n",
    "        xw, w = apply_window(s, window_name)\n",
    "        f, mag = one_sided_fft(xw, fs)\n",
    "        mag = mag / coherent_gain(w)\n",
    "\n",
    "        if fmax is not None:\n",
    "            valid = f <= fmax\n",
    "            f = f[valid]\n",
    "            mag = mag[valid]\n",
    "\n",
    "        rows.append(mag_to_db(mag))\n",
    "        if freq_axis is None:\n",
    "            freq_axis = f\n",
    "\n",
    "    S = np.array(rows).T\n",
    "    time_axis = sts + segment_duration/2\n",
    "    return time_axis, freq_axis, S\n",
    "\n",
    "time_axis, freq_axis, S = spectral_map(\n",
    "    x_noisy, fs,\n",
    "    segment_duration=0.25,\n",
    "    hop_duration=0.05,\n",
    "    window_name=\"hann\",\n",
    "    fmax=400\n",
    ")\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.imshow(\n",
    "    S,\n",
    "    origin=\"lower\",\n",
    "    aspect=\"auto\",\n",
    "    extent=[time_axis[0], time_axis[-1], freq_axis[0], freq_axis[-1]]\n",
    ")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"Frequency (Hz)\")\n",
    "plt.title(\"Manual time-frequency map (dB)\")\n",
    "plt.colorbar(label=\"Magnitude (dB)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "624bdd3f",
   "metadata": {},
   "source": [
    "### Use the time-frequency map to answer:\n",
    "\n",
    "- When does the second tone appear?\n",
    "- At what time does the chirp begin?\n",
    "- Over what frequency range does the chirp evolve?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "613b010d",
   "metadata": {},
   "source": [
    "## 10. Spectrogram\n",
    "\n",
    "This is a compact built-in way to visualize what you just constructed manually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd905ccc",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 5))\n",
    "plt.specgram(x_noisy, NFFT=256, Fs=fs, noverlap=200)\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"Frequency (Hz)\")\n",
    "plt.title(\"Built-in spectrogram\")\n",
    "plt.ylim(0, 400)\n",
    "plt.colorbar(label=\"Intensity (dB)\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.14.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}