分数阶傅里叶变换 (FRFT) 原理与MATLAB实现

分数阶傅里叶变换(FRFT)是传统傅里叶变换的广义形式，在信号处理领域具有重要应用。

一、FRFT数学原理

1.1 基本定义

分数阶傅里叶变换定义为：

\(\mathcal{F}^α[f(t)](u) = \int_{-\infty}^{\infty} K_α(t,u)f(t)dt\)

其中核函数为：

\(K_α(t,u) = \begin{cases} \sqrt{\frac{1-j\cot\alpha}{2\pi}} \exp\left(j\frac{t^2+u^2}{2}\cot\alpha - jtu\csc\alpha\right) & \alpha \neq n\pi \\ \delta(t-u) & \alpha = 2n\pi \\ \delta(t+u) & \alpha = (2n+1)\pi \end{cases}\)

1.2 重要性质

性质	数学表达式
线性性	\(\mathcal{F}^α[af(t)+bg(t)] = a\mathcal{F}^α[f] + b\mathcal{F}^α[g]\)
旋转相加性	\(\mathcal{F}^{α+β} = \mathcal{F}^α\mathcal{F}^β\)
时移特性	\(\mathcal{F}^α[f(t-\tau)] = e^{j\frac{\tau^2}{2}\sin\alpha\cos\alpha} e^{ju\tau\sin\alpha}\mathcal{F}^α[f](u-\tau\cos\alpha)\)
调制特性	\(\mathcal{F}^α[e^{jωt}f(t)] = e^{j\frac{ω^2}{2}\sin\alpha\cos\alpha} e^{-juω\cos\alpha}\mathcal{F}^α[f](u-ω\sin\alpha)\)
能量守恒	\(\int|f(t)|^2dt = \int|\mathcal{F}^α[f](u)|^2du\)

1.3 离散FRFT计算

采用Ozaktas采样型算法，计算复杂度O(N log N)：

F^α = D^{-1}P^αDf

其中：

• \(D\)：傅里叶变换矩阵
• \(P^α\)：分数阶相位矩阵

二、MATLAB实现代码

2.1 核心FRFT函数

function F = frft(f, a)
% 分数阶傅里叶变换
% 输入:
%   f - 输入信号
%   a - 变换阶数 (0 ≤ a ≤ 4)
% 输出:
%   F - 分数阶傅里叶变换结果

% 确保a在[0,4]范围内
a = mod(a, 4);

% 特殊情况处理
if a == 0
    F = f;
    return;
elseif a == 1
    F = fftshift(fft(fftshift(f)));
    return;
elseif a == 2
    F = flipud(f);
    return;
elseif a == 3
    F = fftshift(ifft(fftshift(f)));
    return;
end

% 常规情况计算
N = length(f);
shft = rem((0:N-1) + fix(N/2), N) + 1;
s = f(shft);

% 计算调频率
alpha = a * pi/2;
cot_alpha = cot(alpha);
csc_alpha = csc(alpha);

% 计算相位因子
exp1 = exp(-1i * pi * tan(alpha/2) * (-N:1:N-1).^2 * (N/2)^(-1));
exp2 = exp(1i * pi * csc_alpha * (-N:1:N-1).^2 * (N/2)^(-1));
exp3 = exp(-1i * pi * tan(alpha/2) * (-N:1:N-1).^2 * (N/2)^(-1));

% Ozaktas算法核心步骤
chrp1 = s .* exp1(N+1:2*N);
c = fft(chrp1);
chrp2 = c .* exp2(N+1:2*N);
d = ifft(chrp2);
F = d .* exp3(N+1:2*N);

% 调整幅度
F = F * sqrt(abs(csc_alpha)/N);
F = F(shft);
end

2.2 可视化工具

function plot_frft_results(t, f, F, a)
% 绘制FRFT结果
% 输入:
%   t - 时间向量
%   f - 原始信号
%   F - FRFT结果
%   a - 变换阶数

figure('Position', [100, 100, 1200, 800]);

% 原始信号
subplot(3,2,1);
plot(t, real(f), 'b', 'LineWidth', 1.5);
title('原始信号(实部)');
xlabel('时间 (s)');
ylabel('幅度');
grid on;

subplot(3,2,2);
plot(t, imag(f), 'r', 'LineWidth', 1.5);
title('原始信号(虚部)');
xlabel('时间 (s)');
ylabel('幅度');
grid on;

% FRFT结果
u = linspace(-pi, pi, length(F)); % 分数阶频率域

subplot(3,2,3);
plot(u, abs(F), 'k', 'LineWidth', 1.5);
title(sprintf('FRFT幅度谱 (α=%.2fπ)', a/2));
xlabel('分数阶频率');
ylabel('幅度');
grid on;

subplot(3,2,4);
plot(u, angle(F), 'm', 'LineWidth', 1.5);
title(sprintf('FRFT相位谱 (α=%.2fπ)', a/2));
xlabel('分数阶频率');
ylabel('相位 (rad)');
grid on;

% 时频平面表示
subplot(3,1,3);
imagesc(t, u, repmat(abs(F), length(t), 1));
colormap('hot');
colorbar;
title(sprintf('分数阶傅里叶平面 (α=%.2fπ)', a/2));
xlabel('时间 (s)');
ylabel('分数阶频率');
axis xy;
end

2.3 应用示例：线性调频信号分析

%% FRFT应用示例：线性调频信号分析
clear; clc; close all;

% 参数设置
fs = 1000;          % 采样频率 (Hz)
T = 1;              % 信号时长 (s)
N = fs * T;         % 采样点数
t = linspace(-T/2, T/2, N); % 时间向量

% 生成线性调频信号
f0 = 10;            % 起始频率 (Hz)
f1 = 100;           % 终止频率 (Hz)
k = (f1 - f0)/T;    % 调频率

% 信号模型: s(t) = exp(j2π(f0*t + 0.5*k*t^2))
s = exp(1j*2*pi*(f0*t + 0.5*k*t.^2));

% 计算最佳变换阶数
alpha_opt = 2*atan2(1, k*T^2)/pi;

% 计算不同阶数的FRFT
a_values = [0, 0.5, alpha_opt, 1, 1.5, 2];
results = cell(length(a_values), 1);

for i = 1:length(a_values)
    a = a_values(i);
    F = frft(s, a);
    results{i} = F;
    
    % 绘制结果
    figure('Name', sprintf('α = %.2fπ', a/2), 'NumberTitle', 'off');
    plot_frft_results(t, s, F, a);
end

%% 最佳阶数能量聚焦分析
F_opt = frft(s, alpha_opt);

% 计算能量集中度
[~, max_idx] = max(abs(F_opt));
energy_window = 50; % 窗口大小
energy_total = sum(abs(F_opt).^2);
energy_focused = sum(abs(F_opt(max_idx-energy_window:max_idx+energy_window)).^2);
concentration_ratio = energy_focused / energy_total;

fprintf('最佳变换阶数: α = %.4fπ\n', alpha_opt/2);
fprintf('能量集中度: %.2f%%\n', concentration_ratio*100);

%% 时频分析对比
figure('Position', [100, 100, 1400, 600]);

% STFT分析
subplot(1,2,1);
window = hamming(128);
noverlap = 120;
nfft = 1024;
spectrogram(s, window, noverlap, nfft, fs, 'yaxis');
title('短时傅里叶变换 (STFT)');
colorbar off;

% FRFT分析
subplot(1,2,2);
u = linspace(-pi, pi, length(F_opt));
imagesc(t, u, repmat(abs(F_opt), length(t), 1));
colormap('hot');
title(sprintf('分数阶傅里叶变换 (α=%.4fπ)', alpha_opt/2));
xlabel('时间 (s)');
ylabel('分数阶频率');
colorbar;
axis xy;

三、FRFT特性分析

3.1 阶数对变换结果的影响

变换阶数 α	物理意义	信号表现
0.0	恒等变换	保持原信号
0.5	1/4阶变换	时频中间状态
αopt	最佳阶数	能量高度集中
1.0	傅里叶变换	频域表示
1.5	3/4阶变换	时频中间状态
2.0	时间反转	信号反转

3.2 线性调频信号处理效果

对于线性调频信号：

s(t) = e^{j2\pi(f_0t + \frac{1}{2}kt^2)}

当选择最佳变换阶数：

\alpha_{opt} = \frac{2}{\pi}\tan^{-1}\left(\frac{1}{kT^2}\right)

FRFT能将信号能量高度集中在分数阶频率域的某一点上，实现最优能量聚集。

四、应用场景

4.1 雷达信号处理

%% 应用：雷达多目标检测
% 参数设置
fs = 1e6;           % 采样率 1MHz
T = 0.1;            % 脉冲宽度 100ms
t = -T/2:1/fs:T/2;  % 时间向量

% 生成三个不同参数的LFM信号
s1 = exp(1j*2*pi*(50*t + 500*t.^2));  % 目标1
s2 = exp(1j*2*pi*(100*t + 800*t.^2)); % 目标2
s3 = exp(1j*2*pi*(150*t + 300*t.^2)); % 目标3

% 加性高斯白噪声
SNR = 10; % 信噪比(dB)
signal = s1 + s2 + s3;
signal_noisy = awgn(signal, SNR, 'measured');

% FRFT检测
alpha_range = linspace(0, 2, 500); % 阶数范围
detection_plane = zeros(length(alpha_range), length(signal));

for i = 1:length(alpha_range)
    F = frft(signal_noisy, alpha_range(i));
    detection_plane(i, :) = abs(F);
end

% 检测结果可视化
figure('Position', [100, 100, 1200, 800]);
imagesc(t, alpha_range, detection_plane);
colormap('jet');
colorbar;
xlabel('时间 (s)');
ylabel('变换阶数 α');
title('多目标FRFT检测平面');

% 标记目标位置
hold on;
scatter(0, 2*atan2(1, 500*T^2)/pi, 100, 'wx', 'LineWidth', 2);
scatter(0, 2*atan2(1, 800*T^2)/pi, 100, 'wx', 'LineWidth', 2);
scatter(0, 2*atan2(1, 300*T^2)/pi, 100, 'wx', 'LineWidth', 2);
legend('目标1', '目标2', '目标3');

4.2 通信系统应用

%% 应用：FRFT域多载波调制
% 参数设置
M = 16;             % 16-QAM调制
N_sub = 64;         % 子载波数量
alpha = 1.2;        % 分数阶变换阶数

% 生成随机数据
data = randi([0 M-1], N_sub, 1);
mod_data = qammod(data, M, 'UnitAveragePower', true);

% 传统OFDM调制
ofdm_signal = ifft(mod_data, N_sub);

% FRFT-OFDM调制
frft_signal = frft(mod_data, alpha);

% 信道传输（多径衰落）
h = [0.8, 0, 0, 0, 0.3, 0, 0, 0.2]; % 多径信道
ofdm_rx = conv(ofdm_signal, h, 'same');
frft_rx = conv(frft_signal, h, 'same');

% 接收端处理
% OFDM解调
ofdm_demod = fft(ofdm_rx);
ofdm_data = qamdemod(ofdm_demod, M, 'UnitAveragePower', true);

% FRFT-OFDM解调
frft_demod = frft(frft_rx, -alpha); % 使用负阶数进行逆变换
frft_data = qamdemod(frft_demod, M, 'UnitAveragePower', true);

% 误码率计算
ofdm_ber = sum(data ~= ofdm_data) / N_sub;
frft_ber = sum(data ~= frft_data) / N_sub;

fprintf('传统OFDM误码率: %.4f\n', ofdm_ber);
fprintf('FRFT-OFDM误码率: %.4f\n', frft_ber);

% 信号星座图对比
figure('Position', [100, 100, 1000, 400]);
subplot(1,2,1);
scatter(real(ofdm_demod), imag(ofdm_demod), 50, 'filled');
title(sprintf('传统OFDM星座图 (BER=%.4f)', ofdm_ber));
grid on;
axis square;

subplot(1,2,2);
scatter(real(frft_demod), imag(frft_demod), 50, 'filled');
title(sprintf('FRFT-OFDM星座图 (BER=%.4f)', frft_ber));
grid on;
axis square;

参考代码 分数阶FRFT变换 youwenfan.com/contentcnp/80516.html

五、性能优化技巧

5.1 计算加速方法

function F = frft_fast(f, a)
% 加速版FRFT实现
% 使用矩阵预计算和内存优化

persistent cot_table csc_table; % 持久化变量保存查表

% 参数检查
if nargin < 2
    error('需要两个输入参数');
end

N = length(f);
a = mod(a, 4);

% 特殊角度快速处理
if a == 0
    F = f;
    return;
elseif a == 1
    F = fftshift(fft(fftshift(f)));
    return;
elseif a == 2
    F = flipud(f);
    return;
elseif a == 3
    F = fftshift(ifft(fftshift(f)));
    return;
end

% 创建/更新三角函数查表
if isempty(cot_table) || size(cot_table, 1) ~= N
    cot_table = zeros(N, 1001);
    csc_table = zeros(N, 1001);
    alpha_vals = linspace(0, 4, 1001);
    
    for i = 1:1001
        alpha = alpha_vals(i) * pi/2;
        cot_table(:, i) = cot(alpha);
        csc_table(:, i) = csc(alpha);
    end
end

% 查找最近的alpha索引
alpha_idx = round(a * 250) + 1; % 1001个点覆盖[0,4]范围
if alpha_idx > 1001
    alpha_idx = 1001;
end

cot_alpha = cot_table(1, alpha_idx);
csc_alpha = csc_table(1, alpha_idx);

% 后续计算与标准frft相同...
% [此处省略与标准frft相同的计算步骤]
end

5.2 数值稳定性改进

1. 避免奇异点：

   % 添加小偏移量避免cot(0)和csc(0)问题
   epsilon = 1e-10;
   if abs(sin(alpha)) < epsilon
       alpha = alpha + epsilon;
   end
   

2. 高精度计算：

   % 使用符号计算提高精度
   if abs(a - 1) < 0.01
       syms t u alpha_sym;
       K_sym = sqrt((1-1i*cot(alpha_sym))/(2*pi)) * ...
               exp(1i*(t^2 + u^2)/2*cot(alpha_sym) - 1i*t*u*csc(alpha_sym));
       % 将符号表达式转换为函数句柄
       K_func = matlabFunction(K_sym, 'Vars', [t, u, alpha_sym]);
       % 使用K_func计算FRFT...
   end
   

六、总结

分数阶傅里叶变换作为傅里叶变换的广义形式，在时频分析、信号检测和通信系统中有广泛应用。本文提供的MATLAB实现具有以下特点：

1. 完整实现：包含核心FRFT算法和可视化工具
2. 应用导向：提供雷达检测和通信系统应用实例
3. 性能优化：包含加速算法和数值稳定性改进
4. 全面分析：展示不同变换阶数下的信号特性变化

通过合理选择变换阶数，FRFT能有效解决传统傅里叶变换在处理非平稳信号（如线性调频信号）时的局限性，实现信号能量的最优聚集。