分数阶傅里叶变换在MATLAB中的实现与chirp信号检测

简介：分数傅里叶变换（FRFT）是传统傅里叶变换在角度域上的推广，可对信号进行连续阶次的时-频分析。这份代码包面向信号处理与图像加密方向的研究者、工程师及高年级学生，提供4种典型MATLAB实现，帮助理解FRFT原理并直接用于科研实验。压缩包共9个文件，包含8个.m源码与1个说明文档，整体仅9KB，代码轻量易读，便于按需选用、扩展及教学演示。现已吸引1915人学习下载。实现方式涵盖直接矩阵运算、基于FFT的快速算法、迭代近似法及窗函数插值法，各有适用场景；结合chirp2rot、freq_shear等辅助函数，可完成时域剪切、尺度变换等操作，适合教学演示与非平稳信号分析，也能对比各方法的精度与计算效率。亦可通过调整阶数α实现图像加密与解密模拟，是兼顾理论验证与工程落地的实用参考。

1. 分数阶傅里叶变换不是玄学：一条chirp信号困住了普通FFT

一条线性调频信号，带宽 200MHz，脉宽 10μs，你把它丢进 FFT，频谱摊成一条宽平顶，峰值高度比理论值掉了快两个数量级。这不是 FFT 错了，是线性调频的瞬时频率一直在变，能量在频域里“站不住”。分数阶傅里叶变换（FRFT）做的事情，是把时频平面旋转一个角度，让这条斜线竖起来，能量在某个分数阶域里聚成一个峰。简单说，普通 FFT 是固定旋转 90°，FRFT 允许你旋转任意角度，这个角度就是阶数 p 在控制。这篇博文要讲的是 MATLAB 里怎么落地：最小可用的函数怎么写、chirp 信号怎么验证、阶数怎么扫、归一化怎么调、以及我踩过的那几个坑。适合做雷达信号处理、光学衍射模拟、时频分析或语音调频检测的工程师。

2. 先弄懂FRFT在算什么：旋转角、阶数p与三种离散化路线

2.1 傅里叶变换是特例，FRFT是时频平面上的旋转

分数阶傅里叶变换的定义可以理解为：把信号从时间轴投影到一条与时间轴夹角为 α 的轴上。这里 α = p·π/2，p 就是阶数。p=0 时变换结果就是原信号本身，p=1 时退化为普通傅里叶变换，p=−1 时是逆傅里叶变换，p=2 时信号会时间反转。所以 FRFT 不是“另一种变换”，它是把 FFT 嵌入到一个连续旋转族里。

数学上，连续 FRFT 的核函数写成

K_α(t,u) = sqrt(1 − i·cot α) · exp(i·π·(cot α·(t²+u²) − 2·csc α·t·u))

当 α = π/2 时 cot α = 0、csc α = 1，核退化成 exp(−i·2π·t·u)，这就是 DFT 的核。这也是验证代码正确性的一个关键点：p=1 时你的 FRFT 函数必须和 fft() 的结果一致，如果有出入，说明离散化缩放或相位符号写错了。

工程上我最常问自己的问题是：这个变换结果应该放在哪个轴上看？普通频谱横轴是频率，分数阶域的横轴既不是时间也不是频率，而是一条旋转后的坐标轴。你只需要记住一个直觉： FRFT 的结果是信号在这条旋转轴上的能量分布，旋转角由 p 决定。

2.2 p从0到1之间发生了什么：chirp信号是如何在分数阶域聚焦的

一条线性调频信号在时频平面上是一条斜线，在普通频域里看它是摊开的。换个角度审视：如果旋转角恰好和这条斜线的倾角互补，那么这条斜线在旋转后的坐标轴上就是一个垂直于轴的点，能量集中成一个峰。这就是为什么 chirp 检测是 FRFT 最经典的应用场景。

关键点在于：最优阶数 p 和 chirp 的调频斜率 k、采样时长 T 直接相关。理论上可以推导出公式，但实际工程里我更习惯直接扫描 p。扫描的做法是：对同一段信号，从 p=0 以固定步长跑到 p=2，对每个 p 做一次 FRFT，取输出幅度最大值，画一条“峰值能量 vs 阶数”的曲线，峰值对应的 p 就是最优阶数。整个过程像调收音机旋钮找台，简单粗暴但非常可靠。

需要注意的是 p 的范围可以超出 [0,1]。p=1.3 和 p=−0.7 在结果上有对称关系，因为旋转角 α 超过 π/2 后进入“过旋转”。实际处理 chirp 时最优 p 经常落在负数区，所以扫描范围我一般取 [−0.5, 1.5]，步长 0.01，找到粗值后用小范围抛物线插值细化。

2.3 三种离散化思路，别一上来就选快速算法

FRFT 的连续定义在 MATLAB 里不能直接算，必须离散化。常见路线有直接采样法、快速分解法和特征分解法，三者的取舍差异很大。

直接采样法就是把积分写成求和，核矩阵 N×N 一次性构造出来，然后用一次矩阵乘法完成变换。它慢，但公式和连续定义一一对应，出错了容易查。快速分解法是把 FRFT 拆成两次 chirp 乘法和一次 FFT，但不同文献对采样间隔、缩放系数的处理不一样，拿到手必须用逆变换自检才能确认没改坏。特征分解法构造离散 Hermite 函数作为特征向量，性质最干净，但实现工作量最大。

我的习惯是：第一次跑通概念用直接采样法，确认物理过程对了之后换成快速分解法做性能测试。如果你一上来就抄一段快速算法的代码，结果不对时你会完全不知道该调哪里。

3. 在MATLAB里把FRFT跑起来：最小函数、chirp验证与逆变换自检

3.1 直接采样法实现disfrft：代码短、概念准、适合验证

下面是我一直保留在工具箱里的直接采样实现。它不求快，求的是“和连续定义对得上”。

function X = disfrft(x, p) % disfrft - 直接采样实现的离散分数阶傅里叶变换 % X = disfrft(x, p) % 输入: x - 输入信号，列向量，长度 N % p - FRFT 阶数，常用范围 [-1, 2] % 输出: X - 分数阶傅里叶域结果，长度 N N = length(x); alpha = p * pi / 2; % 两个边界情况直接退化处理 if abs(alpha) < 1e-12 X = x(:); % p=0: 原信号 return; end if abs(alpha - pi/2) < 1e-12 X = fft(x(:)); % p=1: 退化为FFT return; end % 时频平面缩放：让采样间隔乘积为 1/N scale = 1 / sqrt(N); t = (0:N-1)' * scale; % 时间采样点 u = t.'; % 旋转域采样点 % 构造核矩阵 K(k,m) a = cot(alpha); b = csc(alpha); T2 = t.^2 * ones(1, N); U2 = ones(N, 1) * (u.^2); phase = pi * (a * (T2 + U2) - 2 * b * (t * u)); K = sqrt(1 - 1i*a) * exp(1i * phase) * scale; X = K * x(:); end

这段代码里最关键的是 scale = 1/sqrt(N)。它让时间轴和旋转域的采样间隔乘积等于 1/N，这样当 p=1 时核矩阵的相位项退化成 exp(−i·2π·m·k/N)，和 DFT 核一致。如果你把 scale 去掉，p=1 时输出和 fft() 之间会差一个常数因子，后续所有幅度比较都会出错。

参数方面，x 传行向量也能跑，因为函数内部统一做了 x(:)。p 的小数位数决定了你的扫描精度，0.01 步长足够工程用。N 建议先用 256 验证，N=1024 时核矩阵已经占 16MB 内存，N=4096 时占 512MB，这时候直接采样法就不合适了，换快速版本。

注意：这个函数每次调用都会重新构造核矩阵。如果你在循环里扫描 p，性能会很差。验证性使用没问题，批量扫描请拉到 3.3 的快速版本。

3.2 用chirp信号验证“旋转”的含义：扫描阶数找峰值

有了 disfrft 之后，先用一条已知参数的 chirp 验证整个链路。下面这段生成一个线性调频信号，然后扫描阶数找峰值。

% 生成线性调频测试信号 fs = 1000; % 采样率 1000 Hz N = 1024; % 采样点数 t = (0:N-1)' / fs; f0 = 50; % 起始频率 50 Hz k = 200; % 调频斜率 200 Hz/s x = exp(1i * 2*pi * (f0 * t + 0.5 * k * t.^2)); % 扫描阶数 p，从 -0.5 到 1.5，步长 0.01 p_list = -0.5:0.01:1.5; peak = zeros(size(p_list)); for i = 1:length(p_list) Xp = disfrft(x, p_list(i)); peak(i) = max(abs(Xp)); end % 找峰值对应的阶数 [pmax, idx] = max(peak); p_opt = p_list(idx); fprintf('最优阶数 p = %.2f，峰值 = %.2f\n', p_opt, pmax);

这段代码核心是“对每个 p 做一次 FRFT，取输出模的最大值”。普通 FFT 里这个 chirp 的频谱是平的，但 FRFT 在某个 p 附近会出现一个尖锐的峰，峰值高度可能比 FFT 谱高 5～20 倍。如果你扫描出来的曲线只有一个明显的峰，说明你的核矩阵公式、缩放系数、信号参数四条链路全对。

实际工程里有两点要提醒。第一，扫描步长 0.01 足够粗定位，但最优 p 大概率落在两个网格点之间，需要抛物线插值细化，后面 4.2 会给代码。第二，上面这个循环每次都重建核矩阵，N=1024 时跑完大概几十秒到一两分钟，属于正常范围；如果你想验证得更快，先用 N=256 的信号。

3.3 快速版本frft_fast：性能翻倍但有前提

直接采样法的瓶颈是核矩阵 O(N²) 的构造和乘法。工程上我一般用快速分解法替代，思路是把 FRFT 拆成“chirp 乘法 → FFT → chirp 乘法”三段。

function X = frft_fast(x, p) % frft_fast - 基于chirp分解的快速近似FRFT % X = frft_fast(x, p) % 与 disfrft 接口一致，复杂度 O(N log N) N = length(x); alpha = p * pi / 2; if abs(alpha) < 1e-12 X = x(:); return; end if abs(alpha - pi/2) < 1e-12 X = fft(x(:)); return; end n = (0:N-1)'; % 第一步：输入信号乘一个chirp c1 = exp(-1i * pi * tan(alpha/2) * n.^2 / N); y1 = x(:) .* c1; % 第二步：与卷积核做循环卷积（补零到2N用FFT加速） M = 2 * N; m = (0:M-1)'; c2 = exp(1i * pi * csc(alpha) * (m - N).^2 / N); y2 = ifft(fft([y1; zeros(N,1)]) .* fft(c2)); y2 = y2(N:2*N-1); % 截取中心对齐段 % 第三步：再乘同一条chirp，整体缩放 y3 = y2 .* c1; X = sqrt(1 - 1i*cot(alpha)) * y3 / sqrt(N); end

这段代码的正确性依赖两个约定：chirp 相位里除以 N，以及最后乘 sqrt(1−i·cot α)/sqrt(N)。不同文献可能把 /N 放在不同的位置，导致输出幅度差一个常数。我的经验是：跑自检脚本，看逆变换残差，不要凭“看起来像”就认为代码没问题。

用快速版本替换 disfrft 后，扫描 p 的循环速度能提升两个数量级，N=1024 时从几十秒降到零点几秒。代价是你必须无条件信任那三个归一化系数。如果不放心，把 3.2 的验证脚本里的 disfrft 换成 frft_fast，对比同一批 p 的峰值曲线，两条曲线应该重合。

4. 分数阶傅里叶变换的避坑清单：四个必踩的坑与排查方法

4.1 阶数p和旋转角α换算错了，整个结果全是乱的

现象：p=1 时 frft 结果和 fft() 对不上，差得离谱；或者扫描 p 时峰值落在完全不合理的位置，比如 0.5 附近出现第二个莫名其妙的峰。原因：把旋转角写成了 α = p·π，而不是 p·π/2。这个错误在代码里非常隐蔽，因为 p=1 时两种写法都给出 π 或 π/2，差 2 倍，只有扫描 p 不为整数的区间时才会暴露。解决：在函数入口处加一行断言，或者做一个 p=1 的特判分支指向 fft()，像我 3.1 和 3.3 里写的那样。还有一个关联坑：FRFT 对阶数有周期性，p 和 p+4 结果相同，p 和 2−p 有对称关系。我见过有人把 1.6 阶的结果当 0.4 阶用，峰值确实还在，但坐标轴方向反了，后续参数全错。

4.2 采样率和chirp速率不匹配，扫描p永远找不到单峰

现象：扫描 p 得到的峰值曲线不只是单峰，而是在一大段 p 上都有较高的铺底，或者出现多个同高度的峰。原因有两个。一是 chirp 的调频斜率 k 太大，导致瞬时频率在窗口内超出了 fs/2，产生混叠。二是窗口内的周期数太少，信号本质上是“一段截断正弦”，FRFT 无法把它压成一个点。解决：先估算 k·T，要求 k·T 在数值上比 fs/2 小一个量级以上；不行就分帧处理或者降低采样率。还有一个细化技巧：找到粗峰后，用峰附近三点做抛物线插值。

% p_opt 粗值两侧各取一点，抛物线插值细化 i0 = idx; p1 = p_list(i0-1); p2 = p_list(i0); p3 = p_list(i0+1); a1 = peak(i0-1); a2 = peak(i0); a3 = peak(i0+1); h = p2 - p1; delta = h * (a1 - a3) / (2 * (a1 - 2*a2 + a3)); p_fine = p2 + delta; fprintf('细化后的最优阶数 p = %.4f\n', p_fine);

这个插值公式在三点抛物线近似的假设下能给出亚步长精度。实际工程里我把步长设为 0.01，插值后 p 的误差通常在 ±0.005 以内，足够反推 chirp 斜率了。

4.3 逆变换不回来：归一化系数和共轭转置的坑

现象：做完 X = frft(x, p) 后，用 y = frft(X, −p) 逆变换，y 和 x 差一个常数倍数，甚至形状相同但时间顺序反了。原因：直接采样法的核矩阵只是连续 FRFT 核的数值近似，并不是严格酉矩阵；另一个常见错误是逆变换时用了 p 而不是 −p，因为“转置”和“共轭转置”在复指数核下是两个不同操作。解决：先检查你用的 p 符号，正变换 p、逆变换必须 −p；再检查归一化，把输入输出都归一化到单位能量后对比。下面这段自检脚本我每次拿到新函数都会先跑一遍。

% 逆变换自检脚本 x = randn(1024,1) + 1i*randn(1024,1); p = 0.3; X = frft_fast(x, p); y = frft_fast(X, -p); err = norm(y - x) / norm(x); fprintf('相对残差 = %.4e\n', err);

如果残差在 1e-2～1e-3 量级，对这个离散化方案来说算正常；如果到 1e-1，说明缩放系数有问题；如果峰值位置还发生了偏移，多半是核矩阵的相位符号反了。追求严格酉性的场景，用特征分解法实现的 FRFT，那类代码自检残差能到 1e-12。

4.4 边界效应与频谱混叠：输出两端的镜面假峰

现象：FRFT 输出在两端出现对称的假峰，看起来像真实峰值在边界处被“反弹”了。原因：输入信号在窗口两端不连续，FRFT 核的周期性相当于把信号不断延拓，边界跳变处产生泄漏能量；另一个来源是信号带宽超出了离散采样体系的约束，在高阶 p 附近被混叠回低频。解决：给输入加窗，Hann 窗足够，让两端平滑落到 0；其次把 N 取大一些，让信号主体远离边界。还有一个工程习惯：扫描 p 时不要看输出整条线，只看中间 70% 区间的最大值，两侧各 15% 直接丢弃，可以有效规避边界假峰的干扰。这个习惯牺牲了一点检测灵敏度，但换来的是稳定性，批量处理数据时非常值得。

5. 让FRFT真正能用的三个进阶习惯

5.1 拿到任何frft代码，第一件事是跑逆变换残差

不管是自己写的还是从别人那里拿的，先不拿真实信号测试，用噪声信号跑逆变换自检，残差不达标就说明归一化或符号有问题。这个检查 10 秒钟能跑完，但能省下后面两小时的排查时间。我的习惯是 p 取 0.3、0.8、1.2 三个值都跑一遍，因为不同 p 下逼近误差不一样。

5.2 二维FRFT：图像和光学问题怎么扩维

FRFT 具有可分离性，二维变换就是先对每一行做一次一维变换，再对每一列做一次一维变换。

function F2 = frft2(img, p_row, p_col) % 二维FRFT：先按行变换，再按列变换 F2 = frft_fast(img.', p_row).'; F2 = frft_fast(F2, p_col); end

光学衍射传播模拟里，这个变换等价于在分数阶域描述波前传播，p_row 和 p_col 可以不同。图像加密场景也常用二维 FRFT，密钥就是两个阶数。要注意的是 frft_fast 处理矩阵时按列操作，所以首次变换前先转置，变换完再转置回来，否则行列被交换。

5.3 大N时的内存与时间管理：核矩阵的爆内存问题

直接采样法在 N=4096 时核矩阵是 4096×4096 个复数，占 512MB 内存，N=8192 直接超过 2GB，一般的机器会卡死。所以大 N 场景我只用 frft_fast，它的内存占用是 O(N) 而不是 O(N²)。另外，我一般先用 N=256 验证公式和流程，确认物理问题理解正确后再上 N=4096 跑最终结果，中间绝不跳级。这个习惯帮我挡掉了至少三次“看起来在跑，其实在死循环”的尴尬。

我在这件事上最大的教训是：FRFT 的代码包看起来都一样，差一个归一化系数结果就完全不能用。别迷信某个现成脚本，拿到手先跑逆变换残差，再跑 chirp 峰值扫描，两道验证都过了才敢往正式流程里放。希望你也能在分数阶傅里叶变换这条路上少走弯路，希望帮到你。

本文还有配套的精品资源，点击获取