OFDM信号PARP的CCDF图

文章目录

- 引言
- 代码
- 代码疑难解答
- 参考文献

引言

本书主要参考了文献1，但实际上该书中符号和表述的错误非常多（只能说棒子是这样的）；同时因为发表时间的关系，很多MATLAB代码进行了更新，原书提供的代码已经无法正常运行。我将修补后的代码给出，并且给出了疑难问题的理解。

代码

% plot_CCDF.m
% Plot the CCDF curves of Fig. 7.3.

clear all; clc; clf
N_OFDM = 2.^[6:10];  %N_OFDM代表OFDM中做IFFT点的数量
b=2; M=2^b; 
Nblk = 1e4; % Nblk Number of Blocks 
PARP_dB = [4:0.1:10]; 
N_PARP_dB = length(PARP_dB);

CCDF_formula=inline('1-((1-exp(-z.^2/(2*s2))).^N)','N','s2','z'); % Eq.(7.9) 构造一个内联函数对象

for n = 1:length(N_OFDM)    

    N=N_OFDM(n); 
    x = zeros(Nblk,N); 
    sqN=sqrt(N);

    for k = 1:Nblk
       X = mapper(b,N); % 生成 N 个 2^b QAM modulated symbols 
       x(k,:) = ifft(X,N)*sqN; % 确保IFFT变换后能量一致
       CFx(k) = PAPR(x(k,:)); % 计算信号的PAPR
    end

    % 通过锐利信道特性计算平均的σ^2的值
    sigma2 = mean(mean(abs(x)))^2/(pi/2);

    % 计算CCDF的理论值  the maximum of Zn is equivalent to the crest factor, sqrt(PAPR)
    CCDF_theoretical=CCDF_formula(N,sigma2,10.^(PARP_dB/20)); 

    for i = 1:N_PARP_dB
       CCDF_simulated(i) = sum(CFx>PARP_dB(i))/Nblk;
    end

    % 绘制对数图
    semilogy(PARP_dB,CCDF_theoretical,'k-');  hold on; grid on;
    semilogy(PARP_dB(1:3:end),CCDF_simulated(1:3:end),'k:*');
end

axis([PARP_dB([1 end]) 1e-2 1]);  % 确认坐标轴范围
title('OFDM system with N-point FFT');
xlabel('PAPR0[dB]'); ylabel('CCDF=Probability(PAPR>PAPR0)'); 
legend('Theoretical','Simulated');

function [modulated_symbols,Mod] = mapper(b,N)
% If N is given, it generates a block of N random 2^b-PSK/QAM modulated symbols.
% Otherwise, it generates a block of 2^b-PSK/QAM modulated symbols for [0:2^b-1].


M=2^b; % Modulation order or Alphabet (Symbol) size

% 生成一个相移键控 PSK 调制器对象
if b==1, Mod='BPSK'; A=1; mod_object=comm.PSKModulator('ModulationOrder', M);
elseif b==2, Mod='QPSK';  A=1;
    %QPSK 调制的信号星座图是一个单位圆的 4 个点，表示 0 度、90 度、180 度和 270 度的相位
    mod_object=comm.PSKModulator('ModulationOrder', M,'PhaseOffset',pi/4);
else
    % 生成一个 QAM 调制器对象
    Mod=[num2str(2^b) 'QAM']; Es=1; A=sqrt(3/2/(M-1)*Es);
    mod_object=comm.RectangularQAMModulator('ModulationOrder', M, 'SymbolMapping', 'Gray');
end

% 虽然这里是用了A来做功率归一化，实际上调制器函数本身就自带有归一化功能
if nargin==2 % generates a block of N random 2^b-PSK/QAM modulated symbols
    modulated_symbols = A*mod_object(randi([0 M-1], N, 1));
else
    modulated_symbols = A*mod_object([0:M-1]');
end

function [PAPR_dB, AvgP_dB, PeakP_dB] = PAPR(x)
% PAPR_dB  : PAPR[dB]
% AvgP_dB  : Average power[dB]
% PeakP_dB : Maximum power[dB]

%MIMO-OFDM Wireless Communications with MATLAB㈢   Yong Soo Cho, Jaekwon Kim, Won Young Yang and Chung G. Kang
%2010 John Wiley & Sons (Asia) Pte Ltd

Nx=length(x); xI=real(x); xQ=imag(x);
Power = xI.*xI + xQ.*xQ;
PeakP = max(Power); PeakP_dB = 10*log10(PeakP);
AvgP = sum(Power)/Nx; AvgP_dB = 10*log10(AvgP);
PAPR_dB = 10*log10(PeakP/AvgP);

代码疑难解答

为什么要使用 x(k,:) = ifft(X,N)*sqN; 这个IFFT公式

这是为了是变换前后时域和频率域的能量保持一致。

$\sigma^2$ 这个式子是如何计算，这个计算是必须的吗？

根据锐利分布的公式，平均值 $E(X)=\sigma\sqrt{\pi/2}$ 。详细推导可以参考概率分布的教科书。但实际上， $\sigma^2$ 根据默认的设置应该是 $\sigma^2=0.5$ 。因为我们对输入功率在mapper 函数中就进行了归一化处理。sigma2 = mean(mean(abs(x)))^2/(pi/2); 这一行代码实际上是没有必要的。下图是设置 $\sigma^2=0.5$ 后仿真结果图。