项目背景详细介绍
循环卷积（Circular Convolution，有时称为周期卷积）在信号处理、数字滤波、频域分析、通信、多速率系统与某些工程算法（例如基于环缓冲的实时滤波、循环相关、循环卷积定理在FFT实现下的加速）中非常常见。与线性卷积不同，循环卷积假设信号在给定长度 NNN 上是周期延拓的——也就是说，索引超出范围的样本将被模 NNN 映射回区间内。这使得循环卷积可以通过长度为 NNN 的离散傅里叶变换（DFT）非常方便地在频域完成：利用卷积定理，周期卷积在频域等于对应 DFT 值的逐点乘积。
在工程应用中，循环卷积有两类常见使用场景：
-
信号尺寸固定为 NNN，且处理要求为周期卷积（例如环形缓冲、循环相关、循环卷积用于实时系统）。
-
线性卷积通过适当的零填充转化为在更大长度上的循环卷积（例如把两个长度分别为 L 与 M 的序列零填充到 N≥L+M−1，然后用长度 N 的 FFT 做点乘再 IFFT，从而得到线性卷积结果的前 L+M−1 项）。这是 FFT 用于加速线性卷积的基础。
教学上实现并对比朴素循环卷积与基于 FFT 的循环卷积，可以帮助学生和工程师直观理解：
-
DFT 与 IDFT 的实现与细节（特别是迭代位倒置 Cooley–Tukey）；
-
频域乘法如何实现卷积（以及为什么要按长度 NNN 进行零填充或截断）；
-
数值误差（虚部残留、舍入误差）与调试方法（阈值裁剪）；
-
实际工程中的性能权衡（小 NNN 用朴素卷积即可，长序列建议 FFT）。
本项目目标是：用纯 Java（无第三方库）实现教学级别的循环卷积算法，包括朴素循环卷积与基于 FFT 的快速循环卷积实现，配套示例、正确性校验与性能比较，适合博客/课堂演示。代码要易读、注释详细、所有源码放在单一代码块（便于复制），并在文后提供方法作用的简明解读，方便直接作为博客内容或课堂讲义。
项目需求详细介绍
功能性需求：
-
实现 朴素循环卷积（直接定义实现），适用于两个长度都为 N 的序列，时间复杂度 O(N^2)。
-
实现 基于 FFT 的循环卷积（使用 DFT/IDFT 的循环卷积定理），支持任意指定长度 NNN（当需要时应要求 N≥max(lenA,lenB)N \ge \max(lenA, lenB)N≥max(lenA,lenB) 并且常为二的幂以便 FFT 高效实现）。
-
提供将线性卷积通过零填充转为循环卷积的示例（即如何通过选择 N≥L+M−1N \ge L+M-1N≥L+M−1 将线性卷积用 FFT 计算）。
-
实现一个教学用的 Complex 复数类和 FFT 类（迭代 Cooley–Tukey，位倒置重排，正/逆变换）。
-
提供
Main示例：包含小规模快速验证（朴素 vs FFT 比较）、大规模性能粗测、以及说明数值误差与阈值处理。
非功能性需求：
-
纯 Java 标准库实现（
java.lang、java.util、java.math等）。 -
易于扩展（例如后续可以替换为高性能 FFT 库或并行版本）。
-
对异常参数（如 N 非法、数组为空）做合理检测并提示有意义错误。
相关技术详细介绍
-
因此可以通过 FFT（快速傅里叶变换）和 IFFT（逆变换）实现周期卷积。
-
FFT 实现要点
-
常用 Cooley–Tukey 算法：迭代实现通常包含位判置（bit-reversal）与多层蝶形运算（butterfly）。
-
长度 NNN 通常取 2 的幂以简化实现与提高效率（但一般算法也可支持任意甚至素因子分解的长度）。
-
逆变换通常在计算完成后对结果除以 NNN 做归一化。
-
-
线性卷积通过 FFT 的实现
若想得到长度 L+M−1 的线性卷积，将序列 a、b 零填充到长度 N≥L+M−1（最好选择 N 为 2 的幂），再分别做 FFT，逐点相乘，最后做 IFFT，取实部并截取前 L+M−1 项即可得到线性卷积。 -
数值误差处理
-
频域计算使用浮点数（double）会引入舍入误差，逆 FFT 后理论上的虚部应接近零，实际会有小量残余（例如 1e−121e^{-12}1e−12 到 1e−91e^{-9}1e−9）；
-
可设置阈值（例如 1e−91e^{-9}1e−9）将接近零的虚部或实部误差裁剪；
-
对于要求很高精度的场景，可以考虑用高精度库或整数 FFT（NTT）在模数域内计算（受限于模数选择与序列中值范围）。
-
-
性能权衡
-
朴素循环卷积：时间 O(N2)O(N^2)O(N2)，当 NNN 很小（如几十、几百）且频繁变换 N 小时可接受；
-
FFT 实现：时间 O(NlogN)，当 N 较大时（通常几百以上）更有优势，但有常数开销（复数运算、三角函数预计算等）。
-
实现思路详细介绍
设计上分为以下模块（教学友好、便于拆分）：
-
Complex：复数类，提供加减乘除、共轭、标量缩放、模值等基础运算。为了降低 GC 压力，尽量在 FFT 核心循环使用原地更新并复用数组对象。 -
FFT：实现迭代 Cooley–Tukey FFT（含位倒置），提供fft(Complex[] a, boolean invert)。若invert==true执行逆变换并在最后归一化。 -
CircularConvolution：提供三类接口：-
naiveCircularConvolution(double[] a, double[] b)：朴素实现，用于长度相同的序列； -
fftCircularConvolution(double[] a, double[] b, int N)：在给定周期长度 NNN 下使用 FFT 实现循环卷积（如果传入 N 小于数组长度抛出异常）；方法会把输入复制或截断成长度 N； -
linearConvolutionViaFFT(double[] a, double[] b)：示例如何将线性卷积转为周期卷积（选择 N >= L+M-1）；
-
-
Utils：辅助函数（如 nextPowerOfTwo、数组复制/填充、最大绝对差比较）。 -
Main：演示与测试：小 N 验证（朴素与 FFT 结果对比并检查误差）、较大 N 性能对比（计时）、并演示线性卷积通过 FFT 的方法。
实现注意事项：
-
fftCircularConvolution要求输入数组长度 <= N；若其中一个长度小于 N 自动零填充；若不等长行为文档化（建议先明确 N）。 -
FFT内部计算三角函数时可用Math.cos/Math.sin，或预计算旋转因子（wlen）以减少重复计算。 -
输出结果取实部（
C[i].re），并在需要时对虚部做阈值判断。 -
代码在博客中以单块形式给出，分注释区段模拟不同文件，便于读者复制到 IDE 并拆分成真实文件。
完整实现代码
// =======================================================
// File: Complex.java
// Description: 教学用复数类（支持基本运算）
// =======================================================
package com.example.circconv;
public class Complex {
public double re; // 实部
public double im; // 虚部
public Complex() { this.re = 0.0; this.im = 0.0; }
public Complex(double re) { this.re = re; this.im = 0.0; }
public Complex(double re, double im) { this.re = re; this.im = im; }
// 加法（返回新对象）
public Complex add(Complex o) {
return new Complex(this.re + o.re, this.im + o.im);
}
// 减法（返回新对象）
public Complex sub(Complex o) {
return new Complex(this.re - o.re, this.im - o.im);
}
// 乘法（返回新对象）
public Complex mul(Complex o) {
return new Complex(this.re * o.re - this.im * o.im, this.re * o.im + this.im * o.re);
}
// 标量乘
public Complex scale(double k) {
return new Complex(this.re * k, this.im * k);
}
// 共轭
public Complex conj() {
return new Complex(this.re, -this.im);
}
// 绝对值（模）
public double abs() {
return Math.hypot(re, im);
}
@Override
public String toString() {
if (im >= 0) return String.format("%.6f+%.6fi", re, im);
else return String.format("%.6f%.6fi", re, im);
}
}
// =======================================================
// File: FFT.java
// Description: 迭代版 Cooley–Tukey FFT（含逆变换）
// Notes:
// - 输入数组长度必须是 2 的幂；若不是，需要先用 Utils.nextPowerOfTwo 填充
// - invert=true 表示计算逆变换（并在末尾做除以 n 归一化）
// =======================================================
package com.example.circconv;
public class FFT {
/**
* 原地计算 FFT 或逆 FFT
* @param a 输入/输出复数数组，长度 n（应为 2 的幂）
* @param invert 是否为逆变换
*/
public static void fft(Complex[] a, boolean invert) {
int n = a.length;
// 位倒置重排
int j = 0;
for (int i = 1; i < n; i++) {
int bit = n >> 1;
while ((j & bit) != 0) {
j ^= bit;
bit >>= 1;
}
j ^= bit;
if (i < j) {
Complex tmp = a[i];
a[i] = a[j];
a[j] = tmp;
}
}
// 逐步蝶形计算
for (int len = 2; len <= n; len <<= 1) {
double ang = 2 * Math.PI / len * (invert ? -1 : 1);
Complex wlen = new Complex(Math.cos(ang), Math.sin(ang));
for (int i = 0; i < n; i += len) {
Complex w = new Complex(1, 0);
for (int k = 0; k < len / 2; k++) {
Complex u = a[i + k];
Complex v = a[i + k + len / 2].mul(w);
a[i + k] = u.add(v);
a[i + k + len / 2] = u.sub(v);
w = w.mul(wlen);
}
}
}
// 逆变换时做归一化
if (invert) {
for (int i = 0; i < n; i++) {
a[i].re /= n;
a[i].im /= n;
}
}
}
/**
* 从实数组构造长度为 n 的 Complex 数组（不足部分用 0 填充）
*/
public static Complex[] toComplex(double[] real, int n) {
Complex[] a = new Complex[n];
for (int i = 0; i < n; i++) {
double val = (i < real.length) ? real[i] : 0.0;
a[i] = new Complex(val, 0.0);
}
return a;
}
}
// =======================================================
// File: Utils.java
// Description: 工具函数（nextPowerOfTwo、数组比较等）
// =======================================================
package com.example.circconv;
public class Utils {
// 返回不小于 x 的最小 2 的幂
public static int nextPowerOfTwo(int x) {
if (x <= 1) return 1;
int n = 1;
while (n < x) n <<= 1;
return n;
}
// 计算两个实数组的最大绝对差（用于比较朴素与 FFT 结果）
public static double maxAbsDiff(double[] a, double[] b) {
int n = Math.min(a.length, b.length);
double max = 0.0;
for (int i = 0; i < n; i++) {
double d = Math.abs(a[i] - b[i]);
if (d > max) max = d;
}
if (a.length != b.length) {
double[] longer = a.length > b.length ? a : b;
for (int i = n; i < longer.length; i++) {
double d = Math.abs(longer[i]);
if (d > max) max = d;
}
}
return max;
}
}
// =======================================================
// File: CircularConvolution.java
// Description: 提供朴素循环卷积与基于 FFT 的循环卷积实现
// =======================================================
package com.example.circconv;
import java.util.Arrays;
public class CircularConvolution {
/**
* 朴素循环卷积（要求 a.length == b.length == N）
* 定义：(a ⨂ b)[n] = sum_{k=0..N-1} a[k] * b[(n-k) mod N]
* 时间复杂度 O(N^2)
*/
public static double[] naiveCircularConvolution(double[] a, double[] b) {
if (a == null || b == null) throw new IllegalArgumentException("输入数组不能为空");
int n = a.length;
if (n != b.length) throw new IllegalArgumentException("朴素循环卷积要求两个数组长度相等");
double[] out = new double[n];
Arrays.fill(out, 0.0);
for (int i = 0; i < n; i++) {
double sum = 0.0;
for (int k = 0; k < n; k++) {
int idx = i - k;
if (idx < 0) idx += n; // 模运算：(i-k) mod n
sum += a[k] * b[idx];
}
out[i] = sum;
}
return out;
}
/**
* 基于 FFT 的循环卷积：
* - a、b 可为任意长度，但要求 N >= max(a.length, b.length)
* - 若 a 或 b 长度小于 N，会用 0 填充
* - 要求 N 为 2 的幂以便 FFT 实现高效（如果不是，建议先用 Utils.nextPowerOfTwo 扩展）
*
* 过程：
* 1) 用 Complex[] 表示 a 和 b（长度 N）
* 2) fft(A), fft(B)
* 3) C[k] = A[k] * B[k]
* 4) ifft(C) -> 实部即为循环卷积结果
*/
public static double[] fftCircularConvolution(double[] a, double[] b, int N) {
if (a == null || b == null) throw new IllegalArgumentException("输入数组不能为空");
if (N <= 0) throw new IllegalArgumentException("N 必须大于 0");
if (a.length > N || b.length > N) throw new IllegalArgumentException("N 必须不小于输入数组长度");
// 将实数组转为复数组并零填充到长度 N
Complex[] A = FFT.toComplex(a, N);
Complex[] B = FFT.toComplex(b, N);
// 正向 FFT
FFT.fft(A, false);
FFT.fft(B, false);
// 频域逐点相乘
Complex[] C = new Complex[N];
for (int i = 0; i < N; i++) {
C[i] = A[i].mul(B[i]);
}
// 逆 FFT
FFT.fft(C, true);
// 取实部作为结果
double[] out = new double[N];
for (int i = 0; i < N; i++) {
// 虚部理论为 0，数值上可能有非常小的残差，视需要进行阈值裁剪
out[i] = C[i].re;
}
return out;
}
/**
* 将线性卷积转成循环卷积的示例：对长度 L 和 M 的序列，若选择 N >= L+M-1，
* 则在长度 N 的循环卷积中取前 L+M-1 项即为线性卷积结果。
*/
public static double[] linearConvolutionViaFFT(double[] a, double[] b) {
int L = a.length;
int M = b.length;
int resultLen = L + M - 1;
int N = Utils.nextPowerOfTwo(resultLen); // 推荐选 2 的幂
double[] circ = fftCircularConvolution(a, b, N);
// 取前 resultLen 项作为线性卷积结果
return Arrays.copyOf(circ, resultLen);
}
}
// =======================================================
// File: Main.java
// Description: 演示与测试（正确性验证与简单性能比较）
// =======================================================
package com.example.circconv;
import java.util.Random;
import java.util.Arrays;
public class Main {
public static void main(String[] args) {
testSmallNCompare();
testLinearViaFFT();
performanceTest();
}
// 小规模正确性对比（朴素 vs FFT）
private static void testSmallNCompare() {
System.out.println("=== 小规模正确性对比 ===");
int N = 16; // 小 N，便于验证
double[] a = new double[N];
double[] b = new double[N];
Random rand = new Random(42);
for (int i = 0; i < N; i++) {
a[i] = rand.nextDouble() - 0.5;
b[i] = rand.nextDouble() - 0.5;
}
double[] naive = CircularConvolution.naiveCircularConvolution(a, b);
double[] fft = CircularConvolution.fftCircularConvolution(a, b, N);
System.out.println("输入 a: " + Arrays.toString(a));
System.out.println("输入 b: " + Arrays.toString(b));
System.out.println("朴素循环卷积: " + Arrays.toString(naive));
System.out.println("FFT 循环卷积: " + Arrays.toString(fft));
double maxDiff = Utils.maxAbsDiff(naive, fft);
System.out.println(String.format("最大绝对误差: %.6e", maxDiff));
}
// 演示如何用 FFT 计算线性卷积（通过选择 N >= L+M-1）
private static void testLinearViaFFT() {
System.out.println("\n=== 线性卷积 via FFT 示例 ===");
double[] a = {1, 2, 3, 4};
double[] b = {1, 0, -1};
double[] linear = CircularConvolution.linearConvolutionViaFFT(a, b);
System.out.println("a: " + Arrays.toString(a));
System.out.println("b: " + Arrays.toString(b));
System.out.println("线性卷积（FFT 得到）: " + Arrays.toString(linear));
// 直接计算线性卷积供参考
double[] direct = naiveLinearConvolution(a, b);
System.out.println("线性卷积（朴素直接）: " + Arrays.toString(direct));
System.out.println("最大绝对误差: " + Utils.maxAbsDiff(linear, direct));
}
// 简单性能粗测：比较朴素 O(N^2) 与 FFT O(N log N)
private static void performanceTest() {
System.out.println("\n=== 简单性能测试 ===");
int N = 1 << 14; // 16384，比较明显能体现 FFT 优势
double[] a = new double[N];
double[] b = new double[N];
Random rand = new Random(123);
for (int i = 0; i < N; i++) {
a[i] = rand.nextDouble() - 0.5;
b[i] = rand.nextDouble() - 0.5;
}
// 朴素可能非常慢，这里只做一次并限制 N 不要太大
long t0 = System.currentTimeMillis();
// 为了安全，避免在大 N 上运行朴素（注释掉），只演示 FFT 时间
// double[] naive = CircularConvolution.naiveCircularConvolution(a, b);
long t1 = System.currentTimeMillis();
double[] fft = CircularConvolution.fftCircularConvolution(a, b, N);
long t2 = System.currentTimeMillis();
System.out.println(String.format("朴素 (未运行，以免太慢)。FFT 耗时: %d ms (数组长度 %d)", (t2 - t1), N));
// 为 sanity，检查 fft 的长度与预期相符并展示少量值
System.out.println("FFT 卷积前 10 项示例: " + Arrays.toString(Arrays.copyOf(fft, Math.min(10, fft.length))));
}
// 直接线性卷积（朴素）用于验证（小数组）
private static double[] naiveLinearConvolution(double[] a, double[] b) {
int L = a.length, M = b.length;
int outLen = L + M - 1;
double[] out = new double[outLen];
Arrays.fill(out, 0.0);
for (int i = 0; i < L; i++) {
for (int j = 0; j < M; j++) {
out[i + j] += a[i] * b[j];
}
}
return out;
}
}
代码详细解读
-
Complex类-
构造方法：构造复数对象（带不同参数）。
-
add,sub,mul,scale,conj,abs：提供复数的基本运算，供 FFT 与频域运算使用。
-
-
FFT类-
fft(Complex[] a, boolean invert)：对数组a原地计算 FFT（当invert=false）或逆 FFT（当invert=true），使用迭代 Cooley–Tukey 算法并包含位倒置重排；逆变换时对结果做除以 n 的归一化。 -
toComplex(double[] real, int n)：把实数组转为长度为n的Complex[]，不足位置填 0。
-
-
Utils类-
nextPowerOfTwo(int x)：返回不小于x的最小 2 的幂，通常用于选择 FFT 长度。 -
maxAbsDiff(double[] a, double[] b)：计算两个实数组的最大绝对差，用于验证结果一致性。
-
-
CircularConvolution类-
naiveCircularConvolution(double[] a, double[] b)：朴素实现的循环卷积，要求两个数组长度相同，按定义直接累加乘积，时间复杂度 O(N2)O(N^2)O(N2)。 -
fftCircularConvolution(double[] a, double[] b, int N)：基于 FFT 的循环卷积实现：先将a、b零填充到长度N，各自做 FFT，频域点乘，做逆 FFT 并取实部，返回长度为N的周期卷积结果。适合较大N，复杂度约 O(NlogN)O(N \log N)O(NlogN)。 -
linearConvolutionViaFFT(double[] a, double[] b)：示例性方法，将线性卷积问题转为选择N >= L+M-1的周期卷积问题，然后调用fftCircularConvolution并截取前L+M-1项得到线性卷积结果。
-
-
Main类（示例）-
testSmallNCompare()：生成小尺寸随机序列，比较朴素循环卷积与 FFT 循环卷积的输出并报告最大绝对误差，用于正确性验证。 -
testLinearViaFFT()：演示用 FFT 计算线性卷积（通过零填充与选择合适N），并与直接朴素线性卷积做对比。 -
performanceTest()：在较大的 N 下运行 FFT 方法并报告耗时（朴素方法在大 N 下耗时极大，示例中不运行朴素以避免超时），用于初步展示 FFT 的性能优势。 -
naiveLinearConvolution(double[] a, double[] b)：朴素线性卷积工具函数，仅用于小规模验证。
-
项目详细总结
本文完整实现了循环卷积的两种主要实现方式：朴素 O(N2)O(N^2)O(N2) 实现与基于 FFT 的 O(NlogN)O(N \log N)O(NlogN) 实现，并演示如何通过零填充将线性卷积转成周期卷积以便使用 FFT 加速。实现包括用于教学的 Complex 类与迭代 FFT，代码结构清晰、注释详尽，适合用作博客教学材料或课堂示范。
主旨要点：
-
周期卷积（循环卷积）在频域中可以通过 DFT 的逐点乘法高效实现；
-
线性卷积通过合适的零填充转换为在更长周期上的循环卷积；
-
FFT 在较大长度场景下能显著降低时间复杂度，但有常数开销与数值误差问题；
-
教学实现适合理解算法与调试，生产环境可替换为高性能库或并行版本以获取更好性能。
项目常见问题及解答
Q1：为什么 FFT 后结果会有小的虚部？
A1：由于浮点数舍入误差与 FFT 的三角函数计算，逆 FFT 后理论应为实数，但实际会有微小虚部（如 1e−121e^{-12}1e−12 级）。可将其视为数值噪声并在展示时阈值裁剪（例如把绝对值小于 1e−91e^{-9}1e−9 的值设为 0）。
Q2：选择 N 时如何决策？
A2：若做 循环卷积（周期卷积）并需要周期长度为 N，则应直接使用该 N（且 N ≥ 输入长度）。若用 FFT 做 线性卷积，必须选择 N ≥ L+M-1；为了效率通常再上取最近的 2 的幂，例如 Utils.nextPowerOfTwo(L+M-1)。
Q3：为什么 FFT 要用 2 的幂？能否用其它长度？
A3：Cooley–Tukey 形式的实现对 2 的幂实现最简单且最快。但 FFT 也可实现任意长度（例如分解为小素因子），但实现复杂度与代码量增加。教学实现通常选择 2 的幂。
Q4：什么时候使用朴素循环卷积？
A4：当 N 很小（例如几十以内），或实现简单、对性能要求不高时，朴素实现更便捷且常数开销小。在实时系统中，若 N 非常小并且资源受限也可选择朴素或优化的滑动窗口算法。
Q5：如何减少 FFT 的常数开销？
A5：预计算复指数（旋转因子）、使用原地算法复用对象、减少数组分配、以及在可能时使用本地高性能库（如 JTransforms）。并行化（多线程或 GPU）也能减少实际耗时。
扩展方向与性能优化
-
并行化 FFT：使用 Java 的
ForkJoinPool或并行流对大数组分段并行计算蝶形，提高多核利用率。注意线程开销与任务粒度。 -
使用高性能 FFT 库：例如 JTransforms、FFTW（通过 JNI）等可显著提升性能与数值稳定性。
-
内存与对象复用：避免频繁 new Complex 对象，在循环中复用数组以减少 GC 压力。
-
混合卷积策略：对于不同大小的核，选择直接卷积或重叠-保存（overlap-save）/重叠-添加（overlap-add）方法结合 FFT 处理长流式信号。
-
高精度需求：若数值误差不能容忍，可考虑使用多精度（BigDecimal）或整数域的 NTT（数论变换，需选择合适模数）。
-
边界处理与窄带优化：在实际信号处理中常结合窗口函数、滤波器设计和可分离滤波器等进一步优化性能与效果。
扫一扫
245
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
