numpy_learning_notes.py

# %%
# -----------------------------------------------------------------------------
# IMPORTING LIBRARIES
# -----------------------------------------------------------------------------
# These are the core libraries used throughout this notebook for numerical
# computing, signal processing, and visualization.

# frist of all what is te meaning of this symbol? (#%%) ---> this is for Jupyter/VS Code interactive cells for open each one of thee cell in sprate ways.
# How to run? ----> open your terminal and enter each command below <<for install>>

# 1. Fourier Transform functions from SciPy (optimized for scientific computing):
# - fft: Converts a signal from the time domain to the frequency domain.
# - ifft: Inverse Fast Fourier Transform (converts frequency domain back to time domain).
# - fftfreq: Generates the frequency axis (in Hz) that matches the FFT output.
# - for install: pip install scipy
from scipy.fft import fft, ifft, fftfreq

# 2. Matplotlib: for plotting and visualizing signals and images.
# - for install: pip install matplotlib
import matplotlib.pyplot as plt

# 3. IPython display: for cleaner output rendering in Jupyter/VS Code interactive cells.
# - for install: pip install ipython
from IPython.display import display

# 4. Pillow (PIL): for opening, reading, and manipulating image files.
# - for install: pip install pillow
from PIL import Image

# 5. NumPy: the core library for array operations and numerical computation in Python.
# - for install: pip install numpy
import numpy as np

# 6. Peak detection function from SciPy's signal processing module.
# Used later in this notebook to demonstrate how to locate local maxima in a signal.
# - for install: already included with scipy (see the FFT import above)
#   no separate install needed. If you skipped that step: pip install scipy
from scipy.signal import find_peaks

# %%
# =============================================================================
# NUMPY BASICS: ARRAY CREATION
# =============================================================================

# Convert a standard Python list into a 1D NumPy array.
a = np.array([1, 2, 3])
type(a)

# %%
a

# %%
# Create an array of 3 zeros.
b = np.zeros(3)
b

# %%
# Create a longer array of 10 zeros.
c = np.zeros(10)
c

# %%
# Reshape the array into a 2D matrix with 10 rows and 1 column.
c.shape = (10, 1)
c

# %%
# Create an array of 10 ones.
d = np.ones(10)
d

# %%
# Check the data type of the first element.
# NumPy defaults to float64 for these array-creation functions.
type(d[0])

# %%
# np.empty allocates memory for the array WITHOUT initializing the values.
# The numbers you see are just whatever was already in that memory location —
# never assume they are zero.
e = np.empty(3)
e

# %%
# np.linspace generates N evenly spaced numbers over a given interval.
# Here: 5 numbers evenly spaced between 2 and 10 (inclusive).
f = np.linspace(2, 10, 5)
f

# %%
# A simple 1D array with two elements.
g = np.array([10, 2])
g

# %%
# A standard Python list.
a_list = [1, 2, 3, 4, 5, 6, 7]

# Wrapping the list in an extra pair of brackets makes NumPy treat it as
# a 2D array with 1 row and 7 columns (shape: (1, 7)), instead of a flat 1D array.
h = np.array([a_list])
h

# %%
# Passing the list directly (no extra brackets) creates a standard 1D array
# with shape (7,).
cv = np.array(a_list)
cv

# %%
type(h)

# %%
# A "list of lists" — Python's native way of representing 2D data.
b_list = [[9, 8, 7, 6, 5, 4, 3], [1, 2, 3, 4, 5, 6, 7]]

# Wrapping this in an extra pair of brackets bumps it up to a 3D array,
# with shape (1, 2, 7).
i = np.array([b_list])
i

# %%
# Passing the list of lists directly creates a standard 2D matrix
# with shape (2, 7) — 2 rows, 7 columns.
zx = np.array(b_list)
zx

# %%
# Check the shape (dimensions) of the 3D array.
i.shape

# %%
# Set a random seed so the "random" numbers below are reproducible every time
# this script runs.
np.random.seed(0)

# Generate 6 random integers between 0 (inclusive) and 10 (exclusive).
j = np.random.randint(10, size=6)
j

# %%
# Access the first element (index 0).
j[0]

# %%
# Slicing: elements from index 0 up to (but not including) index 2.
j[0:2]

# %%
# Access the last element using a negative index.
j[-1]

# %%
# Slicing with a step of 2: start at index 0, stop before index 5, step by 2.
j[0:5:2]

# %%
# =============================================================================
# IMAGE PROCESSING WITH NUMPY
# =============================================================================
# Images are just large numerical arrays under the hood. This section treats
# an image as a matrix of pixel values to demonstrate how NumPy's array
# operations apply directly to image data.

# 1. Load an image from disk using PIL.
img = Image.open("image/616151.jpg")

# 2. Convert the image object into a NumPy array (a matrix of pixel values).
# From this point on, the image is pure numbers — every pixel's brightness
# and color channel is now an element in a matrix.
img = np.array(img)

# Print the image dimensions: (Height, Width, Color Channels).
# The "3" for color channels represents Red, Green, and Blue (RGB).
print(img.shape)

# %%
# Attempt to render the array as an image using IPython's display function.
# (Note: depending on your environment, this may just print the raw numeric
# matrix instead of rendering a visual image.)
display(img)

# %%
type(img)

# %%
# Proper way to visualize the image: matplotlib's imshow() takes the numeric
# matrix and renders it back into a visible picture.
plt.imshow(img)
# plt.axis("off")  # Uncomment to hide the axis ticks for a cleaner look.
plt.show()

# %%
# Vertical flip (upside down):
# img[::-1] reverses the order of the first axis (rows = the Y/height axis),
# which flips the image top-to-bottom.
plt.imshow(img[::-1])

# %%
# Horizontal flip (mirror):
# img[:, ::-1] keeps the rows in order but reverses the columns
# (the X/width axis), producing a left-right mirror image.
plt.imshow(img[:, ::-1])

# %%
# Cropping:
# Simply slice the rows and columns you want to keep.
# Here: keep rows 380–800 and columns 0–1000.
plt.imshow(img[380:800, 0:1000])

# %%
# Downsampling (reducing resolution):
# img[::2, ::2] takes every 2nd row and every 2nd column, effectively
# halving the resolution in both directions — a fast way to shrink large images.
plt.imshow(img[::2, ::2])

# %%
img
# Broadcasting: applying a math function across the entire array at once.
# np.sin() computes the sine of every single pixel value in the matrix —
# NumPy does this over millions of elements almost instantly, without
# needing an explicit loop.
img_sin = np.sin(img)
img_sin

# %%
# =============================================================================
# STATISTICAL & MATHEMATICAL FUNCTIONS
# =============================================================================

# Sum of every value (pixel) in the matrix.
print(np.sum(img))

# Product of every value in the matrix (can overflow or hit 0 easily on large images).
print(np.prod(img))

# Mean (average) of all values.
print(np.mean(img))

# Standard deviation of all values.
print(np.std(img))

# Variance of all values.
print(np.var(img))

# The smallest value (pixel intensity) in the entire array.
print(np.min(img))

# The largest value (pixel intensity) in the entire array.
print(np.max(img))

# The INDEX (flat position) of the smallest value.
print(np.argmin(img))

# The INDEX (flat position) of the largest value.
print(np.argmax(img))

# %%
# =============================================================================
# MASKING & FILTERING
# =============================================================================
k = np.array([1, 2, 3, 4, 5])

# Applying a comparison operator to an array returns a Boolean array,
# element by element: is each value greater than 3?
# Result: [False, False, False, True, True]
k > 3

# %%
k < 3

# %%
# Filtering: using a Boolean array as an index returns only the elements
# where the condition was True.
# Result: [4, 5]
k[k > 3]

# %%
# np.where is a powerful conditional function with 3 arguments:
# np.where(condition, value_if_true, value_if_false)
# Here: "If a pixel's intensity is greater than 100, set it to 255 (pure white);
# otherwise, set it to 0 (pure black)." This is called thresholding —
# a common technique for turning a grayscale/color image into a high-contrast
# black-and-white mask.
img_masked = np.where(img > 100, 255, 0)

# %%
# Display the thresholded (high-contrast) image.
plt.imshow(img_masked)

# %%
# A second thresholding example with a lower cutoff and a softer "on" value.
img_masked_2 = np.where(img > 50, 200, 0)
plt.imshow(img_masked_2)

# %%
# =============================================================================
# MATRIX OPERATIONS & BROADCASTING
# =============================================================================
a_array = np.array([1, 2, 3, 4, 5])
b_array = np.array([6, 7, 8, 9, 10])

# %%
# Element-wise addition: matching positions are added together.
a_array + b_array

# %%
# Broadcasting: adding a single scalar to an entire array.
# NumPy automatically applies "+30" to every element individually.
a_array + 30

# %%
# Element-wise multiplication: matching positions are multiplied together.
a_array * b_array

# %%
# Scalar multiplication: every element is multiplied by 10.
a_array * 10

# %%
# Dot product (matrix multiplication) using the @ operator (Python 3.5+):
# (1*6) + (2*7) + (3*8) + (4*9) + (5*10) = 130
a_array @ b_array

# %%
# Combining slicing with the .T (transpose) attribute:
# img[:, :, 0] selects every row and column, but only channel index 0 (Red).
# .T then swaps the row/column axes, which rotates and mirrors the resulting
# single-channel image.
plt.imshow(img[:, :, 0].T)

# %%
m = np.array([3, 5, 7, 2, 4])

# np.sort returns the elements in ascending order: [2, 3, 4, 5, 7]
np.sort(m)

# %%
# =============================================================================
# SIGNAL GENERATION & PLOTTING
# =============================================================================

# Two simple 1D arrays representing x and y coordinates.
x = np.array([1, 2, 3, 4])
y = np.array([10, 20, 25, 30])

# A basic line plot.
plt.plot(x, y)
plt.show()

# %%
# A common beginner mistake: trying to plot a "sine wave" using just a
# handful of manually-picked points. Because there are so few points,
# matplotlib just connects them with straight lines — this produces a
# jagged, triangular shape, NOT a smooth continuous sine wave.
q = np.array([-2, -1, 0, 1, 2, 3, 4])
n = np.array([0, 1, 0, -1, 0, 1, -1])

# %%
plt.plot(q, n)
plt.show()

# %%
# The correct way to plot a smooth sine wave:
# Generate MANY closely-spaced points using np.linspace so the curve
# looks continuous instead of jagged.
# Here: 100 points evenly spaced between 0 and 2 seconds.
w = np.linspace(0, 2, 100)

# Define the signal's frequency (in Hz). This is a 2 Hz signal.
v = 2

# Apply the standard sine wave formula: sin(2 * pi * frequency * time)
l = np.sin(2 * np.pi * v * w)

# %%
plt.plot(w, l)
plt.grid(True)
plt.show()

# %%
# Now let's combine two sine waves of different frequencies — for example,
# a 3 Hz signal and a 50 Hz signal — to simulate a "real-world" scenario
# where a clean low-frequency signal is contaminated by high-frequency noise.

# Define the first signal (3 Hz).
new_1 = np.linspace(-2, 2, 2000)
new_f_1 = 3
new_fuction_1 = np.sin(2 * np.pi * new_f_1 * new_1)

# %%
# Define the second signal (50 Hz).
new_2 = np.linspace(-2, 2, 2000)
new_f_2 = 50
new_fuction_2 = np.sin(2 * np.pi * new_f_2 * new_2)

# Combine the two signals by simple addition — this is what a "noisy" signal
# looks like: a real signal plus unwanted noise, added together sample by sample.
nb = new_fuction_1 + new_fuction_2

# %%
# Plot the combined (noisy-looking) signal.
plt.plot(nb)
plt.grid(True)
plt.show()

# %%
# =============================================================================
# LEARNING THE FAST FOURIER TRANSFORM (FFT)
# =============================================================================
# 1. Build a synthetic signal in the time domain.

# Sampling rate: how many data points we collect per second.
sampling_rate = 100

# Sample spacing: the time gap between two consecutive samples.
# This is simply the inverse of the sampling rate (1/100 = 0.01 seconds).
d = 1.0 / sampling_rate

# Build the time axis (from 0 to 1 second, in steps of "d").
t = np.arange(0, 1.0, d)
n = len(t)  # Total number of samples (100 in this case).

# Build a combined signal made of a 5 Hz wave and a smaller 20 Hz wave.
# Signal = 5 Hz wave + (half-strength) 20 Hz wave
signal = np.sin(2 * np.pi * 5 * t) + 0.5 * np.sin(2 * np.pi * 20 * t)

# Plot the raw signal in the time domain.
plt.plot(t, signal)
plt.title("Original Signal in the Time Domain")
plt.xlabel("Time (seconds)")
plt.ylabel("Amplitude")
plt.grid(True)
plt.show()

# %%
# 2. Compute the FFT:
# This analyzes the combined signal above and breaks it down into the
# individual frequencies that make it up. The output is an array of
# complex numbers (magnitude + phase information for each frequency).
fft_result = np.fft.fft(signal)

# %%
# 3. Build the matching frequency axis using np.fft.fftfreq:
# - First argument (n): total number of samples.
# - Second argument (d): time spacing between samples.
# Output: an array of frequencies (in Hz) corresponding to each FFT value.
freqs = np.fft.fftfreq(n, d)

# Note: fftfreq returns frequencies in a specific, non-intuitive order —
# it starts at 0, goes up through the positive frequencies, then jumps
# to the most negative frequency and counts back up to -1.
freqs

# %%
# 4. Reordering the frequencies with np.fft.fftshift (for a standard,
# human-readable plot):
# By default, fftfreq's output is NOT sorted from smallest to largest.
# fftshift rearranges it so negative frequencies come first, zero sits
# in the middle, and positive frequencies come last — this is the order
# a normal plot expects.
freqs_shifted = np.fft.fftshift(freqs)
fft_shifted = np.fft.fftshift(fft_result)

# Take the magnitude (absolute value) of the complex FFT output to get
# the "strength" of each frequency component.
magnitude = np.abs(fft_shifted)

# Plot the full frequency spectrum (both negative and positive frequencies).
plt.plot(freqs_shifted, magnitude)
plt.title("Full Frequency Spectrum (using fftfreq + fftshift)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.grid(True)
plt.show()

# %%
# 5. The practical, engineering approach — plot only the positive frequencies:
# For real-valued signals, the FFT output is always symmetric: the negative
# frequencies are just a mirror image of the positive ones and carry no
# extra information. So in practice, we usually only plot the positive half.

# Build a Boolean mask that selects only frequencies greater than 0.
pos_mask = freqs > 0

pos_freqs = freqs[pos_mask]
pos_magnitude = np.abs(fft_result)[pos_mask]

# Plot the final, clean spectrum — you should see two clear peaks,
# exactly at 5 Hz and 20 Hz, matching the two frequencies we combined earlier.
plt.plot(pos_freqs, pos_magnitude)
plt.title("Frequency Analysis (peaks exactly at 5 Hz and 20 Hz)")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.grid(True)
plt.show()

# %%
# Back to the earlier combined signal (3 Hz + 50 Hz):
rt = fft(nb)

# To convert the FFT output into a proper frequency axis, we need to tell
# NumPy the sample spacing (d) — here d=0.002 seconds, matching the 2000
# points we generated over a 4-second window (-2 to 2).
wwe = np.fft.fftfreq(len(nb), d=0.002)

# %%
plt.plot(wwe, np.abs(rt))
plt.grid(True)
plt.show()

# %%
# As before, fftfreq gives both positive AND negative frequencies (because
# the FFT of a real-valued signal is symmetric). To see only the positive
# half cleanly, we slice the array in half:
nj = len(nb)
freqs = np.fft.fftfreq(nj, d=0.002)

# %%
plt.plot(freqs[:nj // 2], np.abs(rt)[:nj // 2])
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# %%
# =============================================================================
# LEARNING SCIPY.SIGNAL.FIND_PEAKS
# =============================================================================
# find_peaks() scans a 1D array and returns the INDICES of all local maxima
# (points that are higher than both of their immediate neighbors).
# It's one of the most useful tools in signal processing for locating
# repeating events in a signal — heartbeats, sensor spikes, vibration
# cycles, etc. Here we demonstrate it on a simple synthetic signal, not
# on real physiological data (that full example lives in the separate
# ECG project repository).

# 1. Build a synthetic "noisy wave" signal to search for peaks in:
# a clean 1 Hz sine wave with a bit of random noise added on top,
# so the peaks aren't perfectly uniform (closer to a real-world signal).
np.random.seed(1)
demo_time = np.linspace(0, 5, 500)
demo_signal = np.sin(2 * np.pi * 1 * demo_time) + 0.15 * np.random.randn(len(demo_time))

plt.plot(demo_time, demo_signal)
plt.title("Synthetic Noisy Signal (for find_peaks demonstration)")
plt.xlabel("Time (seconds)")
plt.ylabel("Amplitude")
plt.grid(True)
plt.show()

# %%
# 2. Naive peak detection — no thresholds at all:
# Without any constraints, find_peaks will flag EVERY local bump, including
# tiny fluctuations caused by noise. This is rarely what you actually want.
peaks_naive, _ = find_peaks(demo_signal)
print(f"Peaks found with no filtering at all: {len(peaks_naive)}")

# %%
# 3. Intelligent peak detection — using thresholds:
# - height: the minimum amplitude a point must reach to count as a real peak.
#   This filters out small noise bumps that never get "tall" enough.
# - distance: the minimum number of samples that must separate two peaks.
#   This prevents two small noise wiggles right next to each other from
#   both being counted as separate peaks.
# Since our sine wave has a period of 1 second and our sampling gives us
# 100 samples/second, real peaks should be roughly 100 samples apart —
# so we set "distance" comfortably below that (e.g. 60) to still catch every
# genuine cycle, while "height" is set just above the noise floor.
peaks_filtered, properties = find_peaks(
    demo_signal,
    height=0.5,
    distance=60
)
print(f"Peaks found with height + distance filtering: {len(peaks_filtered)}")
print(f"Peak indices: {peaks_filtered}")

# %%
# 4. Visualize both results side by side to see the difference filtering makes.
plt.figure(figsize=(10, 4))

plt.plot(demo_signal, label="Original Signal", color="gray", alpha=0.7)

plt.plot(
    peaks_naive, demo_signal[peaks_naive],
    "o", color="orange", alpha=0.4, label="Unfiltered Peaks (includes noise)"
)

plt.plot(
    peaks_filtered, demo_signal[peaks_filtered],
    "x", color="red", markersize=10, markeredgewidth=3, label="Filtered Peaks (real cycles)"
)

plt.axhline(y=0.5, color="green", linestyle="--", label="Height Threshold (0.5)")

plt.title("find_peaks(): Effect of height and distance filtering")
plt.xlabel("Sample Index")
plt.ylabel("Amplitude")
plt.legend(loc="upper right")
plt.grid(True, linestyle=":", alpha=0.6)
plt.show()

# %%
# Takeaway:
# - Without filtering, find_peaks treats every small wiggle as a "peak" —
#   including ones caused entirely by noise.
# - With sensible height/distance thresholds, it correctly isolates the
#   real, repeating cycles of the underlying signal.
# This exact technique — tuning height and distance — is what makes
# find_peaks() useful for real-world tasks like detecting heartbeats in an
# ECG signal, counting steps in accelerometer data, or finding spikes in
# sensor readings. (See the separate `biomedical_signal_analysis` project
# for that real-world application on actual ECG data.)