Advanced Features

In this chapter, we'll explore NumPy's advanced features, including random number generation, polynomial operations, Fourier transforms, file I/O, and more. These features make NumPy a powerful tool for scientific computing.

Random Number Generation

Basic Random Number Generation

import numpy as np

# Set random seed for reproducible results
np.random.seed(42)

# Basic random number generation
print("Basic random number generation:")
print(f"Single random number: {np.random.random()}")
print(f"Random integer: {np.random.randint(1, 10)}")
print(f"Random array: {np.random.random(5)}")
print(f"Random integer array: {np.random.randint(1, 10, 5)}")
print()

# Random numbers from different distributions
print("Random numbers from different distributions:")
print(f"Normal distribution: {np.random.normal(0, 1, 5)}")
print(f"Uniform distribution: {np.random.uniform(-1, 1, 5)}")
print(f"Exponential distribution: {np.random.exponential(2, 5)}")
print(f"Poisson distribution: {np.random.poisson(3, 5)}")
print()

# Multidimensional random arrays
print("Multidimensional random arrays:")
random_2d = np.random.random((3, 4))
print(random_2d)
print()

# Random selection from array
array = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
print("Random selection from array:")
print(f"Select 1: {np.random.choice(array)}")
print(f"Select 3: {np.random.choice(array, 3)}")
print(f"Select 3 without replacement: {np.random.choice(array, 3, replace=False)}")
print(f"Weighted selection: {np.random.choice(array, 3, p=[0.1]*10)}")

Modern Random Number Generator (Recommended)

import numpy as np

# Use new random number generator
rng = np.random.default_rng(42)

print("Modern random number generator:")
print(f"Random number: {rng.random()}")
print(f"Normal distribution: {rng.normal(0, 1, 5)}")
print(f"Integers: {rng.integers(1, 10, 5)}")
print()

# Random permutation
array = np.arange(10)
print(f"Original array: {array}")
rng.shuffle(array)  # Shuffle in place
print(f"After shuffle: {array}")

# Generate permutation
original = np.arange(10)
permutation = rng.permutation(original)
print(f"Original array: {original}")
print(f"New permutation: {permutation}")
print()

# Random sampling
data = np.arange(100)
sample = rng.choice(data, size=10, replace=False)
print(f"Random sample: {sample}")

Random Number Applications

import numpy as np

# Monte Carlo method to estimate π
def estimate_pi(n_samples):
    rng = np.random.default_rng(42)
    # Generate random points in unit square
    x = rng.uniform(-1, 1, n_samples)
    y = rng.uniform(-1, 1, n_samples)
    
    # Count points inside unit circle
    inside_circle = (x**2 + y**2) <= 1
    pi_estimate = 4 * np.sum(inside_circle) / n_samples
    
    return pi_estimate

print("Monte Carlo estimation of π:")
for n in [1000, 10000, 100000, 1000000]:
    pi_est = estimate_pi(n)
    error = abs(pi_est - np.pi)
    print(f"Samples: {n:7d}, π estimate: {pi_est:.6f}, Error: {error:.6f}")
print(f"Actual π value: {np.pi:.6f}")
print()

# 2D Random Walk
def random_walk_2d(steps):
    rng = np.random.default_rng(42)
    # Random directions: up, down, left, right
    directions = np.array([[0, 1], [0, -1], [-1, 0], [1, 0]])
    
    # Randomly choose directions
    chosen_directions = rng.choice(4, steps)
    moves = directions[chosen_directions]
    
    # Calculate positions
    positions = np.cumsum(moves, axis=0)
    
    return positions

print("2D Random Walk:")
walk = random_walk_2d(20)
print(f"Starting position: (0, 0)")
print(f"Final position: ({walk[-1, 0]}, {walk[-1, 1]})")
print(f"Maximum distance: {np.max(np.sqrt(np.sum(walk**2, axis=1))):.2f}")
print(f"Average distance: {np.mean(np.sqrt(np.sum(walk**2, axis=1))):.2f}")

Polynomial Operations

Polynomial Basics

import numpy as np

# Create polynomial p(x) = 2x³ + 3x² - x + 5
# Coefficients in descending order
coeffs = [2, 3, -1, 5]
print(f"Polynomial coefficients: {coeffs}")
print("Polynomial: 2x³ + 3x² - x + 5")
print()

# Evaluate polynomial
x_values = np.array([0, 1, 2, -1])
results = np.polyval(coeffs, x_values)
print(f"x values: {x_values}")
print(f"p(x): {results}")
print()

# Verify calculations
for x, result in zip(x_values, results):
    manual_calc = 2*x**3 + 3*x**2 - x + 5
    print(f"p({x}) = {result}, Manual calculation: {manual_calc}")
print()

# Polynomial roots
roots = np.roots(coeffs)
print(f"Polynomial roots: {roots}")
print("Verify roots:")
for i, root in enumerate(roots):
    value = np.polyval(coeffs, root)
    print(f"  Root {i+1}: {root:.4f}, p(root) = {value:.2e}")

Polynomial Operations

import numpy as np

# Define two polynomials
# p1(x) = x² + 2x + 1 = (x + 1)²
# p2(x) = x - 1
p1 = [1, 2, 1]
p2 = [1, -1]

print("Polynomial operations:")
print(f"p1(x) = x² + 2x + 1")
print(f"p2(x) = x - 1")
print()

# Polynomial addition
sum_poly = np.polyadd(p1, p2)
print(f"p1 + p2 = {sum_poly}")
print(f"i.e.: x² + 3x")
print()

# Polynomial subtraction
diff_poly = np.polysub(p1, p2)
print(f"p1 - p2 = {diff_poly}")
print(f"i.e.: x² + x + 2")
print()

# Polynomial multiplication
mul_poly = np.polymul(p1, p2)
print(f"p1 * p2 = {mul_poly}")
print(f"i.e.: x³ + x² - x - 1")
print()

# Polynomial division
quotient, remainder = np.polydiv(p1, p2)
print(f"p1 / p2:")
print(f"  Quotient: {quotient}")
print(f"  Remainder: {remainder}")
print(f"  i.e.: (x² + 2x + 1) = (x - 1) * (x + 3) + 4")
print()

# Polynomial derivative
derivative = np.polyder(p1)
print(f"Derivative of p1: {derivative}")
print(f"i.e.: 2x + 2")
print()

# Polynomial integral
integral = np.polyint(p1)
print(f"Integral of p1: {integral}")
print(f"i.e.: x³/3 + x² + x + C")

Polynomial Fitting

import numpy as np

# Generate noisy data
np.random.seed(42)
x_data = np.linspace(0, 10, 20)
# True function: y = 2x² - 3x + 1
y_true = 2 * x_data**2 - 3 * x_data + 1
y_noisy = y_true + np.random.normal(0, 5, len(x_data))

print("Polynomial fitting:")
print(f"Number of data points: {len(x_data)}")
print(f"True function: y = 2x² - 3x + 1")
print()

# Fit with different degrees
for degree in [1, 2, 3, 5]:
    coeffs = np.polyfit(x_data, y_noisy, degree)
    
    # Calculate goodness of fit
    y_fit = np.polyval(coeffs, x_data)
    r_squared = 1 - np.sum((y_noisy - y_fit)**2) / np.sum((y_noisy - np.mean(y_noisy))**2)
    
    print(f"Degree {degree} polynomial fit:")
    print(f"  Coefficients: {coeffs}")
    print(f"  R²: {r_squared:.4f}")
    
    if degree == 2:
        print(f"  Fit function: {coeffs[0]:.2f}x² + {coeffs[1]:.2f}x + {coeffs[2]:.2f}")
        print(f"  Compare with true: 2x² - 3x + 1")
    print()

# Predict new data points
x_new = np.array([2.5, 7.5])
best_fit_coeffs = np.polyfit(x_data, y_noisy, 2)  # Use quadratic
y_pred = np.polyval(best_fit_coeffs, x_new)
y_true_new = 2 * x_new**2 - 3 * x_new + 1

print("Prediction results:")
for x, pred, true in zip(x_new, y_pred, y_true_new):
    print(f"x = {x}: Predicted = {pred:.2f}, True = {true:.2f}, Error = {abs(pred-true):.2f}")

Fourier Transform

Basic Fourier Transform

import numpy as np

# Create signal
t = np.linspace(0, 1, 1000, endpoint=False)
freq1, freq2 = 5, 10  # 5Hz and 10Hz
signal = np.sin(2 * np.pi * freq1 * t) + 0.5 * np.sin(2 * np.pi * freq2 * t)

print("Fourier Transform Analysis:")
print(f"Signal length: {len(signal)}")
print(f"Sampling frequency: {1/(t[1]-t[0]):.0f} Hz")
print(f"Signal contains frequencies: {freq1} Hz and {freq2} Hz")
print()

# Calculate FFT
fft_result = np.fft.fft(signal)
freqs = np.fft.fftfreq(len(signal), t[1] - t[0])

# Calculate amplitude spectrum
amplitude = np.abs(fft_result)
phase = np.angle(fft_result)

# Consider only positive frequencies
positive_freq_idx = freqs > 0
positive_freqs = freqs[positive_freq_idx]
positive_amplitude = amplitude[positive_freq_idx]

# Find main frequency components
peak_indices = np.where(positive_amplitude > 0.1 * np.max(positive_amplitude))[0]
peak_freqs = positive_freqs[peak_indices]
peak_amplitudes = positive_amplitude[peak_indices]

print("Detected frequency components:")
for freq, amp in zip(peak_freqs, peak_amplitudes):
    if freq < 50:  # Only show low frequency components
        print(f"  Frequency: {freq:.1f} Hz, Amplitude: {amp:.1f}")
print()

# Inverse Fourier Transform
reconstructed = np.fft.ifft(fft_result)
reconstruction_error = np.mean(np.abs(signal - np.real(reconstructed)))
print(f"Reconstruction error: {reconstruction_error:.2e}")
print()

# Frequency domain filtering example
# Create low-pass filter (keep frequencies below 15Hz)
cutoff_freq = 15
filtered_fft = fft_result.copy()
filtered_fft[np.abs(freqs) > cutoff_freq] = 0

# Inverse transform to get filtered signal
filtered_signal = np.real(np.fft.ifft(filtered_fft))

print(f"Low-pass filter (cutoff: {cutoff_freq} Hz):")
print(f"Original signal energy: {np.sum(signal**2):.2f}")
print(f"Filtered energy: {np.sum(filtered_signal**2):.2f}")
print(f"Energy retention: {np.sum(filtered_signal**2)/np.sum(signal**2)*100:.1f}%")

2D Fourier Transform

import numpy as np

# Create 2D signal (simulated image)
x = np.linspace(0, 1, 64)
y = np.linspace(0, 1, 64)
X, Y = np.meshgrid(x, y)

# Create image with different frequency components
image = (np.sin(2 * np.pi * 5 * X) * np.cos(2 * np.pi * 3 * Y) + 
         0.5 * np.sin(2 * np.pi * 10 * X) * np.sin(2 * np.pi * 8 * Y))

print("2D Fourier Transform:")
print(f"Image size: {image.shape}")
print(f"Image statistics: mean={np.mean(image):.3f}, std={np.std(image):.3f}")
print()

# 2D FFT
fft_2d = np.fft.fft2(image)
fft_shifted = np.fft.fftshift(fft_2d)  # Shift zero frequency to center

# Calculate amplitude spectrum
amplitude_spectrum = np.abs(fft_shifted)
log_amplitude = np.log(1 + amplitude_spectrum)

print("Frequency domain analysis:")
print(f"Maximum amplitude: {np.max(amplitude_spectrum):.1f}")
print(f"Spectrum center position: ({image.shape[0]//2}, {image.shape[1]//2})")
print()

# Inverse transform verification
reconstructed = np.real(np.fft.ifft2(fft_2d))
reconstruction_error = np.mean(np.abs(image - reconstructed))
print(f"Reconstruction error: {reconstruction_error:.2e}")
print()

# Frequency domain filtering
# Create low-pass filter
center_x, center_y = image.shape[0]//2, image.shape[1]//2
y_idx, x_idx = np.ogrid[:image.shape[0], :image.shape[1]]
distance = np.sqrt((x_idx - center_x)**2 + (y_idx - center_y)**2)

# Gaussian low-pass filter
sigma = 10
gaussian_filter = np.exp(-(distance**2) / (2 * sigma**2))

# Apply filter
filtered_fft = fft_shifted * gaussian_filter
filtered_image = np.real(np.fft.ifft2(np.fft.ifftshift(filtered_fft)))

print("Gaussian low-pass filter:")
print(f"Filter sigma: {sigma}")
print(f"Original image energy: {np.sum(image**2):.2f}")
print(f"Filtered energy: {np.sum(filtered_image**2):.2f}")
print(f"Energy retention: {np.sum(filtered_image**2)/np.sum(image**2)*100:.1f}%")

File Input/Output

Text File Operations

import numpy as np
import tempfile
import os

# Create sample data
data = np.random.random((5, 3))
print("Original data:")
print(data)
print()

# Save to text file
with tempfile.NamedTemporaryFile(mode='w', suffix='.txt', delete=False) as f:
    txt_filename = f.name
    
np.savetxt(txt_filename, data, fmt='%.6f', delimiter=',')
print(f"Data saved to: {txt_filename}")

# Load from text file
loaded_data = np.loadtxt(txt_filename, delimiter=',')
print("Data loaded from file:")
print(loaded_data)
print(f"Data matches: {np.allclose(data, loaded_data)}")
print()

# Save multiple arrays
array1 = np.array([1, 2, 3, 4, 5])
array2 = np.array([10, 20, 30, 40, 50])
header_info = "Array1,Array2"

combined_data = np.column_stack([array1, array2])
with tempfile.NamedTemporaryFile(mode='w', suffix='.csv', delete=False) as f:
    csv_filename = f.name

np.savetxt(csv_filename, combined_data, fmt='%d', delimiter=',', header=header_info)
print(f"CSV file saved to: {csv_filename}")

# Load CSV file
loaded_csv = np.loadtxt(csv_filename, delimiter=',', skiprows=1)
print("Data loaded from CSV:")
print(loaded_csv)
print()

# Cleanup temporary files
os.unlink(txt_filename)
os.unlink(csv_filename)

Binary File Operations

import numpy as np
import tempfile
import os

# Create test data
array_int = np.array([1, 2, 3, 4, 5], dtype=np.int32)
array_float = np.array([1.1, 2.2, 3.3, 4.4, 5.5], dtype=np.float64)
array_2d = np.random.random((3, 4))

print("Binary file operations:")
print(f"Integer array: {array_int}")
print(f"Float array: {array_float}")
print(f"2D array shape: {array_2d.shape}")
print()

# Save single array to binary file
with tempfile.NamedTemporaryFile(suffix='.npy', delete=False) as f:
    npy_filename = f.name

np.save(npy_filename, array_2d)
print(f"2D array saved to: {npy_filename}")

# Load binary file
loaded_2d = np.load(npy_filename)
print(f"Loaded array shape: {loaded_2d.shape}")
print(f"Data matches: {np.allclose(array_2d, loaded_2d)}")
print()

# Save multiple arrays to one file
with tempfile.NamedTemporaryFile(suffix='.npz', delete=False) as f:
    npz_filename = f.name

np.savez(npz_filename, 
         integers=array_int, 
         floats=array_float, 
         matrix=array_2d)
print(f"Multiple arrays saved to: {npz_filename}")

# Load multiple arrays
loaded_arrays = np.load(npz_filename)
print("Loaded arrays:")
for key in loaded_arrays.files:
    print(f"  {key}: shape={loaded_arrays[key].shape}, dtype={loaded_arrays[key].dtype}")
print()

# Compressed save
with tempfile.NamedTemporaryFile(suffix='.npz', delete=False) as f:
    compressed_filename = f.name

np.savez_compressed(compressed_filename,
                   integers=array_int,
                   floats=array_float,
                   matrix=array_2d)

# Compare file sizes
original_size = os.path.getsize(npz_filename)
compressed_size = os.path.getsize(compressed_filename)
compression_ratio = original_size / compressed_size

print(f"File size comparison:")
print(f"  Original file: {original_size} bytes")
print(f"  Compressed file: {compressed_size} bytes")
print(f"  Compression ratio: {compression_ratio:.2f}")
print()

# Cleanup temporary files
os.unlink(npy_filename)
os.unlink(npz_filename)
os.unlink(compressed_filename)

Memory-Mapped Files

import numpy as np
import tempfile
import os

# Create large array (simulating big data)
large_array = np.random.random((1000, 1000))

print("Memory-mapped files:")
print(f"Large array shape: {large_array.shape}")
print(f"Array size: {large_array.nbytes / 1024 / 1024:.2f} MB")
print()

# Save as memory-mapped file
with tempfile.NamedTemporaryFile(suffix='.dat', delete=False) as f:
    memmap_filename = f.name

# Create memory-mapped file
memmap_array = np.memmap(memmap_filename, dtype='float64', mode='w+', shape=(1000, 1000))
memmap_array[:] = large_array[:]
memmap_array.flush()  # Ensure data is written to disk

print(f"Memory-mapped file created: {memmap_filename}")
print(f"File size: {os.path.getsize(memmap_filename) / 1024 / 1024:.2f} MB")
print()

# Open memory-mapped file in read-only mode
readonly_memmap = np.memmap(memmap_filename, dtype='float64', mode='r', shape=(1000, 1000))

# Verify data
print("Data verification:")
print(f"Original array mean: {np.mean(large_array):.6f}")
print(f"Memory-mapped mean: {np.mean(readonly_memmap):.6f}")
print(f"Data matches: {np.allclose(large_array, readonly_memmap)}")
print()

# Partial read (advantage of memory mapping)
print("Partial data reading:")
subset = readonly_memmap[100:200, 100:200]  # Only read small portion
print(f"Subset shape: {subset.shape}")
print(f"Subset mean: {np.mean(subset):.6f}")
print()

# Cleanup
del memmap_array, readonly_memmap  # Close memory maps
os.unlink(memmap_filename)

Performance Optimization and Memory Management

Memory Usage Optimization

import numpy as np
import sys

print("Memory usage optimization:")

# Impact of data types on memory
size = 1000000

array_float64 = np.ones(size, dtype=np.float64)
array_float32 = np.ones(size, dtype=np.float32)
array_int32 = np.ones(size, dtype=np.int32)
array_int16 = np.ones(size, dtype=np.int16)
array_int8 = np.ones(size, dtype=np.int8)

print(f"Array size: {size:,} elements")
print(f"float64: {array_float64.nbytes / 1024 / 1024:.2f} MB")
print(f"float32: {array_float32.nbytes / 1024 / 1024:.2f} MB")
print(f"int32:   {array_int32.nbytes / 1024 / 1024:.2f} MB")
print(f"int16:   {array_int16.nbytes / 1024 / 1024:.2f} MB")
print(f"int8:    {array_int8.nbytes / 1024 / 1024:.2f} MB")
print()

# View vs Copy
original = np.arange(1000000)
print(f"Original array memory: {original.nbytes / 1024 / 1024:.2f} MB")

# Create view (no extra memory)
view = original[::2]  # Every other element
print(f"View shares memory: {np.shares_memory(original, view)}")
print(f"View base object: {view.base is original}")

# Create copy (extra memory)
copy = original[::2].copy()
print(f"Copy shares memory: {np.shares_memory(original, copy)}")
print(f"Copy base object: {copy.base is None}")
print()

# In-place operations vs creating new arrays
test_array = np.random.random(100000)
original_id = id(test_array)

# In-place operation
test_array += 1
print(f"Same ID after in-place operation: {id(test_array) == original_id}")

# Creating new array operation
test_array = test_array + 1
print(f"Same ID after new array operation: {id(test_array) == original_id}")

Performance Testing and Optimization

import numpy as np
import time

def time_function(func, *args, **kwargs):
    """Measure function execution time"""
    start_time = time.time()
    result = func(*args, **kwargs)
    end_time = time.time()
    return result, end_time - start_time

print("Performance testing and optimization:")

# Create test data
size = 1000000
array1 = np.random.random(size)
array2 = np.random.random(size)

print(f"Test array size: {size:,}")
print()

# Test different calculation methods
print("Calculate sqrt(x² + y²) using different methods:")

# Method 1: Using np.sqrt
result1, time1 = time_function(lambda: np.sqrt(array1**2 + array2**2))
print(f"Method 1 (np.sqrt): {time1:.4f} seconds")

# Method 2: Using np.hypot
result2, time2 = time_function(lambda: np.hypot(array1, array2))
print(f"Method 2 (np.hypot): {time2:.4f} seconds")

# Method 3: Using np.linalg.norm
stacked = np.stack([array1, array2], axis=0)
result3, time3 = time_function(lambda: np.linalg.norm(stacked, axis=0))
print(f"Method 3 (np.linalg.norm): {time3:.4f} seconds")

print(f"Results equal: {np.allclose(result1, result2) and np.allclose(result2, result3)}")
print(f"Fastest method: Method {np.argmin([time1, time2, time3]) + 1}")
print()

# Memory layout impact on performance
print("Memory layout impact on performance:")

# C-style (row-major) vs Fortran-style (column-major)
array_c = np.random.random((1000, 1000))
array_f = np.asfortranarray(array_c)

print(f"C-style array: {array_c.flags['C_CONTIGUOUS']}")
print(f"Fortran-style array: {array_f.flags['F_CONTIGUOUS']}")

# Row sum
_, time_c_row = time_function(lambda: np.sum(array_c, axis=1))
_, time_f_row = time_function(lambda: np.sum(array_f, axis=1))

# Column sum
_, time_c_col = time_function(lambda: np.sum(array_c, axis=0))
_, time_f_col = time_function(lambda: np.sum(array_f, axis=0))

print(f"C-style row sum: {time_c_row:.4f} seconds")
print(f"C-style column sum: {time_c_col:.4f} seconds")
print(f"F-style row sum: {time_f_row:.4f} seconds")
print(f"F-style column sum: {time_f_col:.4f} seconds")

Practical Examples

Example 1: Digital Signal Processing

import numpy as np

def design_fir_filter(cutoff_freq, sampling_freq, num_taps):
    """Design FIR low-pass filter"""
    # Normalized cutoff frequency
    normalized_cutoff = cutoff_freq / (sampling_freq / 2)
    
    # Generate ideal low-pass filter impulse response
    n = np.arange(num_taps)
    h = np.sinc(2 * normalized_cutoff * (n - (num_taps - 1) / 2))
    
    # Apply Hamming window
    window = np.hamming(num_taps)
    h = h * window
    
    # Normalize
    h = h / np.sum(h)
    
    return h

def apply_filter(signal, filter_coeffs):
    """Apply FIR filter"""
    return np.convolve(signal, filter_coeffs, mode='same')

# Generate test signal
fs = 1000  # Sampling frequency
t = np.linspace(0, 1, fs, endpoint=False)

# Combined signal: 50Hz + 150Hz + noise
signal = (np.sin(2 * np.pi * 50 * t) + 
         0.5 * np.sin(2 * np.pi * 150 * t) + 
         0.2 * np.random.normal(0, 1, len(t)))

print("Digital Signal Processing Example:")
print(f"Sampling frequency: {fs} Hz")
print(f"Signal length: {len(signal)} samples")
print(f"Signal contains: 50Hz + 150Hz + noise")
print()

# Design low-pass filter (cutoff 100Hz)
cutoff = 100
num_taps = 101
filter_coeffs = design_fir_filter(cutoff, fs, num_taps)

print(f"FIR filter design:")
print(f"  Cutoff frequency: {cutoff} Hz")
print(f"  Filter order: {num_taps}")
print(f"  Filter coefficient range: [{np.min(filter_coeffs):.6f}, {np.max(filter_coeffs):.6f}]")
print()

# Apply filter
filtered_signal = apply_filter(signal, filter_coeffs)

# Analyze filtering effect
original_power = np.mean(signal**2)
filtered_power = np.mean(filtered_signal**2)
print(f"Filtering effect:")
print(f"  Original signal power: {original_power:.4f}")
print(f"  Filtered power: {filtered_power:.4f}")
print(f"  Power retention: {filtered_power/original_power*100:.1f}%")
print()

# Frequency domain analysis
original_fft = np.fft.fft(signal)
filtered_fft = np.fft.fft(filtered_signal)
freqs = np.fft.fftfreq(len(signal), 1/fs)

# Calculate attenuation at 50Hz and 150Hz
freq_50_idx = np.argmin(np.abs(freqs - 50))
freq_150_idx = np.argmin(np.abs(freqs - 150))

attenuation_50 = 20 * np.log10(np.abs(filtered_fft[freq_50_idx]) / np.abs(original_fft[freq_50_idx]))
attenuation_150 = 20 * np.log10(np.abs(filtered_fft[freq_150_idx]) / np.abs(original_fft[freq_150_idx]))

print(f"Frequency attenuation:")
print(f"  Attenuation at 50Hz: {attenuation_50:.1f} dB")
print(f"  Attenuation at 150Hz: {attenuation_150:.1f} dB")

Example 2: Advanced Image Processing

import numpy as np

def create_test_image(size=128):
    """Create test image"""
    x = np.linspace(-2, 2, size)
    y = np.linspace(-2, 2, size)
    X, Y = np.meshgrid(x, y)
    
    # Create complex pattern
    image = (np.sin(3 * X) * np.cos(3 * Y) + 
            0.5 * np.sin(8 * X) * np.sin(8 * Y) + 
            0.3 * np.random.random((size, size)))
    
    # Normalize to [0, 255]
    image = ((image - np.min(image)) / (np.max(image) - np.min(image)) * 255).astype(np.uint8)
    return image

def gaussian_kernel(size, sigma):
    """Generate Gaussian kernel"""
    kernel = np.zeros((size, size))
    center = size // 2
    
    for i in range(size):
        for j in range(size):
            x, y = i - center, j - center
            kernel[i, j] = np.exp(-(x**2 + y**2) / (2 * sigma**2))
    
    return kernel / np.sum(kernel)

def convolve_2d(image, kernel):
    """2D convolution"""
    pad_h, pad_w = kernel.shape[0] // 2, kernel.shape[1] // 2
    padded = np.pad(image, ((pad_h, pad_h), (pad_w, pad_w)), mode='reflect')
    
    result = np.zeros_like(image, dtype=np.float32)
    
    for i in range(image.shape[0]):
        for j in range(image.shape[1]):
            result[i, j] = np.sum(padded[i:i+kernel.shape[0], j:j+kernel.shape[1]] * kernel)
    
    return result.astype(np.uint8)

def edge_detection(image):
    """Edge detection"""
    # Sobel operators
    sobel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    sobel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
    
    # Calculate gradients
    grad_x = convolve_2d(image, sobel_x)
    grad_y = convolve_2d(image, sobel_y)
    
    # Calculate gradient magnitude
    magnitude = np.sqrt(grad_x.astype(np.float32)**2 + grad_y.astype(np.float32)**2)
    magnitude = (magnitude / np.max(magnitude) * 255).astype(np.uint8)
    
    return magnitude

# Image processing pipeline
print("Advanced Image Processing Example:")

# Create test image
original_image = create_test_image(64)
print(f"Original image size: {original_image.shape}")
print(f"Pixel value range: [{np.min(original_image)}, {np.max(original_image)}]")
print()

# Gaussian blur
gaussian_5x5 = gaussian_kernel(5, 1.0)
blurred_image = convolve_2d(original_image, gaussian_5x5)

print(f"Gaussian blur:")
print(f"  Kernel size: 5x5")
print(f"  Sigma: 1.0")
print(f"  Blurred pixel range: [{np.min(blurred_image)}, {np.max(blurred_image)}]")
print()

# Edge detection
edges = edge_detection(original_image)
print(f"Edge detection:")
print(f"  Edge pixel range: [{np.min(edges)}, {np.max(edges)}]")
print(f"  Edge pixel count: {np.sum(edges > 50)}")
print()

# Image statistics
print(f"Image statistics comparison:")
print(f"  Original image mean: {np.mean(original_image):.1f}")
print(f"  Blurred image mean: {np.mean(blurred_image):.1f}")
print(f"  Edge image mean: {np.mean(edges):.1f}")
print()
print(f"  Original image std: {np.std(original_image):.1f}")
print(f"  Blurred image std: {np.std(blurred_image):.1f}")
print(f"  Edge image std: {np.std(edges):.1f}")

Chapter Summary

In this chapter, we explored NumPy's advanced features:

• Random Number Generation: Various distributions, random sampling, Monte Carlo methods
• Polynomial Operations: Polynomial calculation, operations, fitting
• Fourier Transform: 1D and 2D FFT, frequency domain analysis and filtering
• File I/O: Text files, binary files, memory-mapped files
• Performance Optimization: Memory management, cache-friendly access, vectorization
• Practical Applications: Signal processing, image processing, and other domain applications

Next Steps

In the next chapter, we'll learn NumPy's practical tips and best practices, including debugging techniques, performance analysis, integration with other libraries, and more.

Exercises

1. Use Monte Carlo method to estimate the definite integral ∫₀¹ x² dx
2. Implement a simple digital filter design function
3. Write an image sharpening filter
4. Use FFT to implement fast convolution
5. Create a memory-efficient large data processing function
6. Implement a polynomial regression class
7. Write a random walk visualization program