从信号处理到机器学习：用Python和NumPy手把手理解傅里叶变换与梯度下降

从信号处理到机器学习：用Python和NumPy手把手理解傅里叶变换与梯度下降

数学理论与工程实践的鸿沟常常让学习者感到困惑——公式推导看似严谨，却难以转化为实际代码。本文将用Python和NumPy作为桥梁，带您体验数学概念从理论到落地的完整过程。我们会发现，傅里叶变换不仅是信号处理的利器，更是理解数据频域特征的窗口；梯度下降也不仅是优化算法的核心，其本质体现了多维空间中的"下山"直觉。

1. 傅里叶变换：时域与频域的双向翻译官

在数字信号处理领域，傅里叶变换就像一位精通两种语言的翻译专家。它能将随时间变化的信号（时域）转换为由不同频率组成的频谱（频域），这种转换保留了原始信号的全部信息，只是换了一种表达方式。

让我们用NumPy生成一个简单的复合信号开始实验：

import numpy as np import matplotlib.pyplot as plt # 生成时间序列 t = np.linspace(0, 1, 1000, endpoint=False) # 创建包含10Hz和20Hz成分的信号 signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 20 * t) plt.plot(t, signal) plt.title('原始时域信号') plt.xlabel('时间(秒)') plt.ylabel('振幅') plt.show()

执行这段代码，你会看到一个波形图，其中包含两个不同频率的正弦波叠加而成。但仅凭时域图像，我们很难准确判断信号中包含哪些频率成分——这正是傅里叶变换大显身手的地方。

1.1 快速傅里叶变换(FFT)实战

NumPy提供了高效的fft函数实现：

# 计算FFT fft_result = np.fft.fft(signal) # 计算频率轴 freqs = np.fft.fftfreq(len(t), t[1]-t[0]) # 绘制频谱图 plt.stem(freqs[:50], np.abs(fft_result)[:50]) plt.title('频域表示') plt.xlabel('频率(Hz)') plt.ylabel('幅度') plt.show()

你会看到两个明显的峰值，分别位于10Hz和20Hz处——这正是我们合成信号时使用的频率。20Hz成分的幅度较小，对应代码中的0.5系数。

注意：实际应用中我们通常只关心正频率部分，且需要进行适当的幅度缩放。完整的处理流程还应包括窗函数选择、零填充等技巧。

1.2 频域滤波的威力

理解了频域表示后，我们可以实现简单的滤波器。比如，要移除信号中的高频噪声：

# 创建带噪声的信号 noisy_signal = signal + 0.2 * np.random.normal(size=len(t)) # 设计低通滤波器 fft_noisy = np.fft.fft(noisy_signal) fft_filtered = fft_noisy.copy() fft_filtered[np.abs(freqs) > 15] = 0 # 截断15Hz以上的频率 # 逆变换回时域 filtered_signal = np.fft.ifft(fft_filtered) plt.figure(figsize=(12,4)) plt.subplot(121) plt.plot(t, noisy_signal) plt.title('带噪声信号') plt.subplot(122) plt.plot(t, filtered_signal.real) # 取实部 plt.title('滤波后信号') plt.show()

这个简单的例子展示了频域操作的强大之处——在时域难以处理的噪声问题，转换到频域后可能变得非常简单。

2. 梯度下降：优化问题的通用解法

从信号处理转向机器学习，梯度下降算法扮演着核心角色。这个看似简单的概念——沿着梯度反方向小步前进——却是训练神经网络等复杂模型的基石。

2.1 梯度下降的几何直观

考虑一个简单的二次函数：

def quadratic(x): return x**2 + 5*x + 6 x = np.linspace(-10, 5, 100) plt.plot(x, quadratic(x)) plt.xlabel('x') plt.ylabel('f(x)') plt.title('二次函数优化') plt.grid() plt.show()

这个函数的导数为f'(x)=2x+5。梯度下降的更新规则为：

def gradient_descent(start, learning_rate, iterations): x = start history = [x] for _ in range(iterations): grad = 2 * x + 5 # 手动计算的导数 x = x - learning_rate * grad history.append(x) return x, history x_opt, history = gradient_descent(start=0, learning_rate=0.1, iterations=20) print(f"找到的最小值点：{x_opt:.4f}") # 可视化优化过程 x_vals = np.array(history) plt.plot(x, quadratic(x)) plt.plot(x_vals, quadratic(x_vals), 'ro-') plt.title('梯度下降轨迹') plt.show()

你会看到算法如何一步步"滑向"函数的最低点。学习率的选择至关重要：

2.2 多元函数的梯度下降

现实中的机器学习问题通常涉及多维参数空间。考虑函数f(x,y)=x²+y²：

def quadratic_2d(x, y): return x**2 + y**2 # 计算梯度 def grad_2d(x, y): return np.array([2*x, 2*y]) # 二维梯度下降 def gd_2d(start, lr, steps): point = np.array(start) path = [point.copy()] for _ in range(steps): grad = grad_2d(*point) point -= lr * grad path.append(point.copy()) return np.array(path) path = gd_2d([5, 5], 0.1, 20) # 绘制等高线和优化路径 x = y = np.linspace(-6, 6, 100) X, Y = np.meshgrid(x, y) Z = quadratic_2d(X, Y) plt.contour(X, Y, Z, levels=20) plt.plot(path[:,0], path[:,1], 'ro-') plt.title('二维梯度下降') plt.show()

这个可视化清晰地展示了算法如何在二维空间中寻找最小值。在更高维度（如神经网络可能有数百万参数）时，原理完全相同，只是难以可视化。

3. 从数学公式到NumPy实现的关键技巧

将数学理论转化为高效代码需要一些技巧。以下是常见数学运算的NumPy实现对照：

3.1 广播机制的应用

NumPy的广播机制能极大简化代码。例如计算L2正则化项：

# 传统实现 def l2_regularization(weights): total = 0.0 for w in weights: total += w**2 return total # NumPy风格 def l2_regularization_np(weights): return np.sum(weights**2)

广播机制同样适用于多维运算：

# 计算一组点到原点的距离 points = np.random.randn(100, 2) # 100个二维点 distances = np.sqrt(np.sum(points**2, axis=1))

3.2 避免显式循环

对于离散傅里叶变换(DFT)：

数学定义： X_k = ∑_{n=0}^{N-1} x_n e^{-i2πkn/N}

低效实现：

def slow_dft(x): N = len(x) X = np.zeros(N, dtype=complex) for k in range(N): for n in range(N): X[k] += x[n] * np.exp(-2j * np.pi * k * n / N) return X

高效向量化实现：

def vectorized_dft(x): N = len(x) n = np.arange(N) k = n.reshape((N, 1)) M = np.exp(-2j * np.pi * k * n / N) return np.dot(M, x)

比较两者的速度差异：

x = np.random.random(1024) %timeit slow_dft(x) # 在我的机器上约4.5秒 %timeit vectorized_dft(x) # 约0.01秒

向量化实现通常能带来数百倍的性能提升，这正是NumPy的强大之处。

4. 机器学习中的数学实践

4.1 线性回归的两种视角

从数学角度看，线性回归可以表示为最小化目标函数：

L(w) = ||Xw - y||²₂

解法一：解析解（正规方程）

w = np.linalg.inv(X.T @ X) @ X.T @ y

解法二：梯度下降

def linear_regression_gd(X, y, lr=0.01, epochs=1000): n_samples, n_features = X.shape w = np.zeros(n_features) for _ in range(epochs): grad = (2/n_samples) * X.T @ (X @ w - y) w -= lr * grad return w

两种方法各有优劣：

4.2 神经网络中的反向传播

反向传播本质上是链式法则的高效实现。考虑一个简单的两层网络：

def relu(x): return np.maximum(0, x) def forward(X, W1, b1, W2, b2): z1 = X @ W1 + b1 a1 = relu(z1) z2 = a1 @ W2 + b2 return z2 def backward(X, y, W1, b1, W2, b2, lr): # 前向传播 z1 = X @ W1 + b1 a1 = relu(z1) z2 = a1 @ W2 + b2 # 反向传播 m = X.shape[0] dz2 = (z2 - y) / m dW2 = a1.T @ dz2 db2 = np.sum(dz2, axis=0) da1 = dz2 @ W2.T dz1 = da1 * (z1 > 0) # ReLU导数 dW1 = X.T @ dz1 db1 = np.sum(dz1, axis=0) # 参数更新 W1 -= lr * dW1 b1 -= lr * db1 W2 -= lr * dW2 b2 -= lr * db2 return W1, b1, W2, b2

这个实现虽然简单，却包含了深度学习的核心数学原理。实际应用中我们会使用优化器（如Adam）和批处理等技巧，但基本框架不变。

5. 性能优化与数值稳定性

5.1 对数域计算

当处理概率或小数值时，直接相乘可能导致下溢：

# 不稳定的实现 def unstable_probability(probs): return np.prod(probs) # 稳定的对数域实现 def stable_log_prob(log_probs): return np.exp(np.sum(log_probs))

5.2 数值稳定的softmax

标准的softmax实现：

def softmax(x): exps = np.exp(x) return exps / np.sum(exps)

数值稳定版本：

def stable_softmax(x): x = x - np.max(x) # 减去最大值防止指数爆炸 exps = np.exp(x) return exps / np.sum(exps)

比较两者的差异：

x = np.array([1000, 1001, 1002]) print(softmax(x)) # 可能得到[nan, nan, nan] print(stable_softmax(x)) # 正确结果

5.3 内存高效的矩阵运算

大矩阵运算时，注意内存布局：

# 低效：创建临时数组 result = np.dot(A, np.dot(B, C)) # 高效：使用einsum result = np.einsum('ij,jk,kl->il', A, B, C)

对于特别大的矩阵，可以考虑分块计算：

def block_matrix_multiply(A, B, block_size=1024): m, n = A.shape _, p = B.shape C = np.zeros((m, p)) for i in range(0, m, block_size): for j in range(0, p, block_size): for k in range(0, n, block_size): ii, jj, kk = slice(i, i+block_size), slice(j, j+block_size), slice(k, k+block_size) C[ii, jj] += A[ii, kk] @ B[kk, jj] return C

在实际项目中，这些优化技巧可能带来显著的性能提升，特别是处理大规模数据时。