OFDM 系统中 PAPR 降低
从零搭一个离散 OFDM 发射链路 → 算 PAPR → 三种"非搜索"降低手段（限幅/压扩/预留子载波）
一、PAPR 的定义（时域）
对一个 OFDM 符号（N 个子载波，sampling pts = N·oversample 的过采样版本）：
PAPR（线性尺度）
工程上常用 dB：
二、主干：OFDM 发射机 + PAPR 计算
2.1 主脚本（ofdm_papr_main.m）
%% ofdm_papr_main.m
% OFDM PAPR 计算 & 降低（不限幅/压扩/预留子载波）
% 不依赖任何搜索/优化工具箱
clear; clc; close all;
rng(0); % 可重复
%% ===== 1. 参数 =====
Nfft = 128; % FFT size（子载波数）
Ncp = 16; % 循环前缀长度
osf = 4; % 过采样倍数（测PAPR需≥4）
Nos = Nfft * osf; % 过采样点数（时域插值法）
Nsymb = 200; % 统计用的OFDM符号数（越多CCDF越稳）
Es = 1; % 平均符号能量归一化（QPSK）
% 子载波分配
% 最简单：直流（k=0）空，两边各用 Nused/2 个数据子载波
Nused = 52; % 实际数据子载波（示意）
k_data = [(1:Nused/2)+1 , Nfft-(Nused/2)+1: Nfft]; % 避开 dc=1和Nyquist
k_reserve = []; % 预留子载波索引（后面TR法用）
fprintf('=== OFDM PAPR 仿真 ===\n');
fprintf('Nfft=%d Ncp=%d osf=%d Nsymb=%d\n',Nfft,Ncp,Nos,Nsymb);
%% ===== 2. 生成一批OFDM符号（原始） =====
tx_syms = cell(Nsymb,1);
time_sigs_orig = zeros(Nos, Nsymb);
for m = 1:Nsymb
% QPSK符号（±1±j）/√2
bits = randi([0 1], Nused, 2);
sym = (1-2*bits(:,1))/sqrt(2) + 1j*(1-2*bits(:,2))/sqrt(2);
X = zeros(Nfft,1);
X(k_data) = sym;
% IFFT → 加CP → 过采样（补零插值 + 低通理想插值用 FFT过采样）
x_base = ifft(X)*sqrt(Nfft); % Nfft点
x_cp = [x_base(end-Ncp+1:end); x_base]; % (Nfft+Ncp)×1
% 过采样：插0再低通 → 用fft技巧
x_os = oversample_by_ifft(x_base, Nfft, osf);
time_sigs_orig(:,m) = x_os;
end
%% ===== 3. PAPR & CCDF =====
[PAPR_orig, ccdf_orig] = papr_ccdf(time_sigs_orig, 0.5);
fprintf('\n原始信号：PAPR0(dB) ~ %.2f dB\n', max(PAPR_orig));
fprintf(' P(PAPR>7dB) ~ %.4f\n', interp1(ccdf_orig.db, ccdf_orig.prob, 7));
%% ===== 4. 方法一：限幅（Clipping）=====
clip_gamma_dB = 6.5; % 限幅门限（相对 RMS），可调
clip_gamma = 10^(clip_gamma_dB/20) * sqrt(mean(abs(time_sigs_orig(:)).^2));
time_clip = cellfun(@(x) clip_signal(x, clip_gamma), ...
num2cell(time_sigs_orig,1), 'UniformOutput',false);
time_clip = cell2mat(time_clip.');
[PAPR_clip, ccdf_clip] = papr_ccdf(time_clip, 0.5);
%% ===== 5. 方法二：μ-law 压扩（Companding）=====
mu = 10; % μ参数，越大压得越狠
time_comp = zeros(size(time_sigs_orig));
for m = 1:Nsymb
x = time_sigs_orig(:,m);
r = abs(x);
r_normed = r / sqrt(mean(r.^2)+eps);
% μ-law 压缩
r_comp = (log(1+mu*r_normed)/log(1+mu)) .* exp(1j*angle(x));
% 恢复平均功率（可选：保证输出功率和输入一致）
gain_corr = sqrt(mean(r.^2)+eps) / sqrt(mean(abs(r_comp).^2)+eps);
time_comp(:,m) = r_comp * gain_corr;
end
[PAPR_comp, ccdf_comp] = papr_ccdf(time_comp, 0.5);
%% ===== 6. 方法三：预留子载波 TR（Tone Reservation）=====
% 思想：选一些子载波永远不送数据，只用来"反相"抵消时域峰值
% 不需要搜索：固定一组随机相位 / 或者解最小二乘（线性系统）→ 也不是搜索
Ntr = 8; % 预留子载波个数
k_tr = randperm(Nfft-2, Ntr) + 1; % 随机挑（避开dc），非搜索
k_tr(k_tr==1)=[]; if numel(k_tr)<Ntr, k_tr(end+1)=Nfft/2; end
k_tr = k_tr(1:Ntr);
time_tr = zeros(Nos, Nsymb);
for m = 1:Nsymb
x_base = time_sigs_orig(1:osf:end, m); % 关键：用基带Nfft点反演
x_base = x_base(1:Nfft); % 确保长度
% 迭代抵消峰值（线性，非搜索）
C = 1.5; iterTR = 5;
for it = 1:iterTR
peaks = find_peaks_for_cancel(x_base, C);
if isempty(peaks), break; end
% 构造 b = -x[peak] 的目标，A = F(:,k_tr) 的列为候选子载波时域基
A = zeros(Nfft, Ntr);
for t = 1:Ntr
A(:,t) = sqrt(Nfft)*ifft(unit_impulse_at(Nfft, k_tr(t)));
end
% 最小二乘（解线性方程，不是搜索）
c = A \ (-x_base(peaks));
% 更新
x_base(peaks) = x_base(peaks) + A(peaks,:)*c;
% 裁剪极端（保底）
x_base = min(abs(x_base),C)*exp(1j*angle(x_base));
end
time_tr(:,m) = oversample_by_ifft(x_base, Nfft, osf);
end
[PAPR_tr, ccdf_tr] = papr_ccdf(time_tr, 0.5);
%% ===== 7. 画图 =====
figure('Color','w','Position',[120 80 1100 420]);
subplot(1,3,1)
plot(ccdf_orig.db, ccdf_orig.prob,'k-','LineW',1.8); hold on
plot(ccdf_clip.db, ccdf_clip.prob,'r--','LineW',1.5)
plot(ccdf_comp.db, ccdf_comp.prob,'b-.','LineW',1.5)
plot(ccdf_tr.db, ccdf_tr.prob, 'g:','LineW',1.8)
grid on; set(gca,'YScale','log')
xlabel('PAPR_0(dB)'); ylabel('CCDF = Prob(PAPR>γ)')
title('CCDF 对比')
legend('原始','限幅','μ压扩','TR预留子载波','Location','southwest')
subplot(1,3,2)
semilogy(ccdf_orig.db, ccdf_orig.prob,'k'); grid on; hold on
semilogy(ccdf_clip.db,ccdf_clip.prob,'r--')
xline(7,'k:'); yline(1e-2,'k:')
title('CCDF (log)')
xlabel('PAPR_0(dB)'); ylabel('Prob')
subplot(1,3,3)
% 抽一个符号看时域包络
x0 = time_sigs_orig(:,1); xc = time_clip(:,1); xt = time_tr(:,1);
plot(abs(x0),'k','LineW',1.2); hold on
plot(abs(xc),'r--','LineW',1.1)
plot(abs(xt),'g:','LineW',1.1)
plot([1 numel(x0)],[clip_gamma clip_gamma],'r:')
grid on; xlabel('样点 n'); ylabel('|x[n]|')
title('一个OFDM符号时域包络对比')
legend('原始','限幅','TR','ClipThr','Location','northwest')
sgtitle('OFDM PAPR 降低方法对比（无搜索工具箱）')
三、底层函数
3.1 过采样（IFFT补零插值）
function x_os = oversample_by_ifft(x_base, Nfft, osf)
% 时域过采样：把Nfft点x放到Nos点IFFT的低位 → 再做一次小IFFT
Nos = Nfft*osf;
X = fft(x_base); % Nfft点频域
X_up = zeros(Nos,1);
% 把正负频搬正确（假设 x_base 是 Nfft 点实对称结构）
% 最稳妥：直流+正频放前面，Nyquist放中间
mid = floor(Nfft/2)+1;
X_up(1:mid) = X(1:mid);
X_up(end-(Nfft-mid)+1:end) = X(mid+1:end);
x_os = ifft(X_up)*sqrt(Nos); % 注意幅度看你是否想保留sqrt(Nfft/Nos)因子
end
3.2 PAPR 计算 + CCDF
function [PAPR_vec, ccdf] = papr_ccdf(time_sigs, rms_ref)
% time_sigs : Nos × Nsymb
% 输出：PAPR_vec每个符号，ccdf结构(db, prob)
Ns = size(time_sigs,2);
PAPR_vec = zeros(Ns,1);
for m = 1:Ns
x = time_sigs(:,m);
p_inst = abs(x).^2;
Pavg = mean(p_inst);
PAPR_vec(m) = max(p_inst)/Pavg;
end
db = 0:0.1:12;
prob = zeros(size(db));
for k = 1:numel(db)
prob(k) = mean(PAPR_vec > 10^(db(k)/10));
end
ccdf.db = db;
ccdf.prob = prob;
end
3.3 限幅（硬限幅，最朴素）
function y = clip_signal(x, gamma)
% gamma = 限幅门限（幅度值，不是dB）
y = x;
mag = abs(x);
over = mag > gamma;
y(over) = gamma .* exp(1j*angle(x(over)));
end
3.4 TR法用到的两个小辅助
function idx = find_peaks_for_cancel(x, C)
% 简单找超过门限C的索引（非搜索优化，只是阈值检测）
idx = find(abs(x) > C);
end
function vec = unit_impulse_at(N,k)
% 频域冲激：第k个bin=1（k从1开始，MATLAB索引）
vec = zeros(N,1);
vec(k) = 1;
end
参考代码 ofdm码中papr的降低 www.youwenfan.com/contentcsv/101627.html
四、重点解释
| 方法 | 核心运算 | 为什么不是搜索 |
|---|---|---|
| 限幅 Clipping | 逐点阈值判断 |x|>γ → 改幅度 |
纯代数/逻辑 |
| μ-law 压扩 | 逐点非线性函数 y=f(|x|,∠x) |
确定性映射 |
| TR 预留子载波 | 每次检测到超门限位置后，解 Ax=b（最小二乘） 更新预留子载波权值 | 解线性方程组（\），不是搜索 |
