初识OFDM(八)：OFDM中的PAPR计算和通频带仿真

初识OFDM(八)：OFDM中的PAPR计算和通频带仿真

零.代码地址

https://github.com/liu-zongxi/OFDM_simulation

代码参考了https://zhuanlan.zhihu.com/p/385096476陈老湿的仿真，但各个函数我都有重新实现，希望写的更规范一些

一. OFDM信号CF的CCFD

1. 代码展示

%-----------------------OFDM信号的PAPR的CCFD----------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年3月31日-----------------%

%%
clear; clc;
Nffts = 2.^[6:10];          % IFFT点数，这是本仿真的变量
Npsk = 2;                   % 调制对应比特数
M= 2^Npsk;                  % 调制阶数
Nblk = 1e3;                 % ?
Z_dBs = 4:0.1:10;           % IFFT后时域信号的能量大小
N_Z_dBs = length(Z_dBs);    % ?
CCDF_formula = @(N, sigma2, z) 1-((1-exp(-z.^2/(2*sigma2))).^N);    % CCFD计算公式

%% 主函数
for kk = 1:length(Nffts)      % 对于每一个不同的Nfft
    Nfft = Nffts(kk);
    x = zeros(Nfft, Nblk);
    x_CF = zeros(1, Nblk);
    CCDF_simulated = zeros(1,N_Z_dBs);
    for ii = 1:Nblk
        X_mod = ModSymbolGenerator(Npsk, Nfft);
        x(:, ii) = ifft(ifftshift(X_mod), Nfft) .*sqrt(Nfft);
        x_CF(ii) = PAPR_dB(x(:,ii));
    end
    sigma2 = mean(mean(abs(x)))^2/(pi/2);   % E(Z^2) = 2sigma^2
    CCDF_theoretical = CCDF_formula(Nfft, sigma2, 10.^(Z_dBs/20));
    for jj = 1:N_Z_dBs
        CCDF_simulated(jj) = sum(x_CF>Z_dBs(jj))/Nblk;
    end
    semilogy(Z_dBs,CCDF_theoretical,'-');  hold on; grid on;
    semilogy(Z_dBs,CCDF_simulated,':*');
end 
title('OFDM system with N-point FFT');
xlabel('PAPR0[dB]');
ylabel('CCDF=Probability(PAPR>PAPR0)'); 
legend('Theoretical','Simulated');

2.代码分析

代码有这么几点是值得注意的

为什么ifft后×sqrt(Nfft)?

x(:, ii) = ifft(ifftshift(X_mod), Nfft) .*sqrt(Nfft);

这是为了保证信号的能量不变，之前已经提过这个问题了

$\sigma$是如何计算的？

这是瑞利分布的特性

由此可以求出${\sigma }^2$

不同变量的dB

我们注意到代码

    CCDF_theoretical = CCDF_formula(Nfft, sigma2, 10.^(Z_dBs/20));

也就是说Z的这个dB表示实际上是$20lo{g}_{10}(Z)$,而能量的dB一般是$10lo{g}_{10}(P)$，这是啥原因呢？不知道有没有大佬能给我解答一下

二.OFDM信号的PDF

1.代码展示

%--------------------------OFDM信道的PDF--------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年3月31日-----------------%

%% 参数设置
clear; clc;
N = 32;          % 生成的调制符号的个数
Npsk = 2;       % 调制的比特
M = 2^Npsk;     % 阶数
Nos = 16;       % 过采样倍数
Nfft = N*Nos;   % fft采样点
Ts = 1/Nfft;    % 采样周期
ts = 0:Ts:1-Ts; % 采样时间
Nhist = 1e3;    % 统计重复次数
Nbin = 30;      % 画直方图时被分成多少类
%% 主程序
for jj = 1:Nhist
    X_mod = ModSymbolGenerator(Npsk, N);
    X(1) = 0+1j*0;
    xI = zeros(Nfft, N);
    xQ = zeros(Nfft, N);
    for ii = 1:N
        % 这是在做什么？计算一下就知道啦
        % 为了作图，每个FFT只使用一个子载波
        % 用那些呢？最左边的和最右边的
        % 为什么？虚拟子载波吧？不太清楚
       if ii<=N/2  
           x = ifft([zeros(1,ii-1) X_mod(ii) zeros(1,Nfft-ii+1)],Nfft);
       else  
           x = ifft([zeros(1,Nfft-N+ii-1) X_mod(ii) zeros(1,N-ii)],Nfft);
       end
       xI(:, ii) = real(x); 
       xQ(:, ii) = imag(x);
    end
    HistI((Nfft*(jj-1)+1):(Nfft*jj)) = sum(xI,2);
    HistQ((Nfft*(jj-1)+1):(Nfft*jj)) = sum(xQ,2);
end
sum_xI = sum(xI,2); sum_xQ = sum(xQ,2); % 这组合出来的就是OFDM符号啦
figure(1), clf, subplot(311);
plot(ts,xI,'k:','linewidth',1), hold on, plot(ts,sum_xI,'b','linewidth',2);
subplot(312)
plot(ts,xQ,'k:','linewidth',1); hold on, plot(ts,sum_xQ,'b','linewidth',2);
subplot(313), plot(ts,abs(sum_xI+1j*sum_xQ),'b','linewidth',2); hold on;
figure(2), clf, subplot(311)
sqrt(var(HistI))
h_xI = histogram(HistI,Nbin);
subplot(312)
h_xQ = histogram(HistQ,Nbin);
subplot(313)
sqrt(var(HistQ))
sqrt(1/(2*N))/Nos
h_x = histogram(abs(HistI+1j*HistI),Nbin);

2.代码分析

这其实是一个很基础的仿真，不过做了这么久的OFDM，怎么能连他的基本的概率特性都不知道呢！所以这个还是很有意义的

如何展现OFDM符号的合成过程

	for ii = 1:N
        % 这是在做什么？计算一下就知道啦
        % 为了作图，每个FFT只使用一个子载波
        % 用那些呢？最左边的和最右边的
        % 为什么？虚拟子载波吧？不太清楚
       if ii<=N/2  
           x = ifft([zeros(1,ii-1) X_mod(ii) zeros(1,Nfft-ii+1)],Nfft);
       else  
           x = ifft([zeros(1,Nfft-N+ii-1) X_mod(ii) zeros(1,N-ii)],Nfft);
       end
       xI(:, ii) = real(x); 
       xQ(:, ii) = imag(x);
    end

这段代码非常适合放在我的第一篇帖子里，他详细的展示了OFDM信号是如何合成的。

此外，为什么要使用两侧的子载波呢？

这也是之前讲过的一个问题了，也就是ifftshift

不过我现在有了一些新的理解

这个循环位移左移动$N/2$，指的是数据的本身长度，而不是补0后的长度，想一想，按这个思路移动后的结果是不是就是中间插0的结果

这个仿真来源于书，他没有前后倒置，显然是错误的，这部分具体可以看我后面的实现

OFDM信号具体是个啥分布啊

这似乎是一个老生常谈的问题，OFDM信号的实部和虚部分别是高斯分布，幅度则是一个瑞利分布

不过，他高斯分布的具体参数是啥呢？这个问题似乎没人说

首先，中心$\mu$肯定是0啦，因为他是cos或者是sin函数的中心极限定理带来的高斯分布

那么方差呢，我数学很差，不过我会仿真啊！结果是$\sigma =\sqrt{\frac{1}{2N}}/Nos$

这个咋理解呢，这实际上是能量开根号啊！然后实部虚部各占一半，不过过采样这个Nos是如何影响的，我目前还没有想清楚

三. 量化的OFDM信号

1.代码展示

%-----------------------OFDM中SNQR的计算-------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年3月31日-----------------%

%%
clear, clf, clc;
Nfft = 64;
Nk = 64;                % 这三个参数不重复说明啦
Npsk = 6;
M = 2^Npsk;
MaxIter=1000;           % 迭代次数
Quantbits = [6:9];      % 量化比特数
V_cutoffs = [2:0.2:8];  % 限幅电平
gss=['ko-';'ks-';'k^-';'kd-'];
SQNRdB = zeros(length(Quantbits), length(V_cutoffs));
%%
for ii = 1:length(Quantbits)
    Quantbit = Quantbits(ii);
    Quantbit_fractional = Quantbit -1;
    for jj = 1:length(V_cutoffs)
        V_cutoff = V_cutoffs(jj);
        Tx = 0; Te = 0;
        for kk = 1:MaxIter
            X_mod = ModSymbolGenerator(Npsk, Nk);
            x = ifft(X_mod, Nfft);
            x = x/sqrt(1/(2*Nk))/V_cutoff;          % 实部或者虚部的方差是(1/(2*Nk))
            % x(:) = (real(x(:))>V_cutoff).*V_cutoff;
%             for kkk= 1:length(x)
%                 if real(x(kkk))>V_cutoff
%                     x(kkk) = V_cutoff + 1j*imag(x(kkk));
%                 end
%                 if imag(x(kkk))>V_cutoff
%                     x(kkk) = real(x(kkk)) + 1j*V_cutoff;
%                 end
%             end
            xq=fi(x,1,Quantbit,Quantbit_fractional);
            xq=double(xq); Px = sum(abs(x).^2); e = x-xq; Pe = sum(abs(e).^2);
            Tx = Tx + Px;  Te = Te + Pe;
        end
        SQNRdB(ii,jj) = 10*log10(Tx/Te)
    end
end
%%
[SQNRdBmax,imax] = max(SQNRdB'); % To find the maximum elements in each row of SQNRdB
for i=1:size(gss,1), plot(V_cutoffs,SQNRdB(i,:),gss(i,:)), hold on;  end
for i=1:size(gss,1)
   str(i,:)=[num2str(Quantbits(i)) 'bit quantization'];
   plot(V_cutoffs(imax(i)),SQNRdBmax(i),gss(i,1:2),'markerfacecolor','r')
end
xlabel('mus(clipping level normalized to \sigma)');  ylabel('SQNR[dB]');
legend(str(1,:),str(2,:),str(3,:),str(4,:));  grid on

2. 代码分析

其实已经没啥分析的了,现在回头看看这个仿真意义不大

不过需要注意的是，他这里的限幅居然不是所谓的cutoff，而是放缩！这真的卡了我很久

四.OFDM的通频带仿真

1.代码展示

%--------------------------OFDM通频带仿真--------------------%
%-----------------------author:lzx-------------------------%
%-----------------------date:2022年4月3日------------------%
%% 设置参数
clear; clf; clc;
% 信噪比
SNRdBs=[0:10];
N_SNR=length(SNRdBs);
% OFDM相关
Maxiter = 1000;          % 信噪比迭代次数
CRs=[0.8:0.2:1.6];      % 限幅系数，他会被用于RMS的相乘限幅
N_CR=length(CRs);
gss='*^<sd>v';
Nmod = 2;               % 映射
M = 2^Nmod;             % 调制阶数
Nk = 128;               % 子载波次数，一般Nk=Nfft
Nfft = 128; 
NGI = 32;               % GI长度
Nframe = 6;             % 一帧6个符号
Nsym = Nfft + NGI;      % 系统长度
% 采样
fs = 1e6;               % 采样次数，也就是两个点直接距离
Nos = 8;                % 过采样倍数
Tsym = 1/(fs/Nsym);     % Nsym占了多长时间？
Ts = 1/(fs*Nos);        % 采样周期
ts = [0:Ts:Nframe*Tsym-Ts].';
% 载波
fc = 2e6;               % 上变频载波频率
wc = 2*pi*fc;           % 角频率
% 滤波器
Fs=8; Norder=104; dens=20; 
FF=[0 1.4 1.5 2.5 2.6 Fs/2];
WW=[10 1 10];
h = firpm(Norder,FF/(Fs/2),[0 0 1 1 0 0],WW,{dens});
% 初始化一些储存数组
CF = zeros(1, Maxiter);
% 一些位置
bers = zeros(N_CR, N_SNR);
% Z相关
Z_dBs = [2:0.1:16];
N_Z = length(Z_dBs);
CCDF = zeros(1, N_Z);
CCDF_CR = zeros(N_CR, N_Z);
%% 主函数
for ii = 1:N_SNR
    SNR = SNRdBs(ii);
    for kk = 1:N_CR
        CR = CRs(kk);
        nobe = 0;
        for jj = 1:Maxiter
            % 生成一帧数据，串并转换，并QPSK，生成一帧
            Frame_FDserial = rand(1,Nk*Nframe*Nmod) > 0.5;     % 发送的是bit
            Frame_FDparallel = reshape(Frame_FDserial,Nk,Nframe*Nmod);% 串并转换
            Frame_mod = QPSKMod(Frame_FDparallel,Nk,Nframe);     %调制
            % 做过采样ifft
            % https://blog.csdn.net/sinat_38151275/article/details/85268026
            Frame_oversampling = [zeros(1, Nframe); Frame_mod(Nk/2+1:Nk,:); zeros(Nk*(Nos-1)-1,Nframe); Frame_mod(1:Nk/2,:)];
            frame = ifft(Frame_oversampling, Nfft*Nos);
%             Frame_mod(1,:) = 0+1j*0;                                    % 没有直流分量
%             frame = IFFTOversampling(Frame_mod, Nk, Nos);           % 过采样IFFT
            frame_GI = AddGI(frame, Nfft*Nos, NGI*Nos, Nframe, "CP");  % 添加保护间隔
            frame_serial = reshape(frame_GI, Nframe*Nsym*Nos, 1);   % 程序写的有点问题，只能都按列来了
            % frame_padding = [zeros((Nfft/2-NGI)*Nos, Nframe); frame_GI; zeros(Nfft*Nos/2, Nframe)];
            frame_passband = sqrt(2) .* real(frame_serial.*exp(1j*wc*ts));
            % frame_clipped = frame_passband;
            frame_clipped = Clipping(frame_passband, CR);
            % frame_filter = frame_clipped;
            frame_filter = ifft(abs(fft(h.',size(frame_clipped, 1))).*fft(frame_clipped));      % 滤波原理
            if ii == N_SNR
                CF(1,jj) = PAPR_dB(frame_filter);
            end
            % 加噪声
            % frame_noise = frame_filter;
            frame_noise = awgn(frame_filter,SNR,'measured');   % add Noise(AWGN)
            frame_baseband = sqrt(2) .* frame_noise.*exp(-1j*wc*ts);
            % frame_sample = frame_baseband(1+[0:Nos:Nsym*Nos*Nframe-1],1);
            frame_parallel = reshape(frame_baseband, Nsym*Nos, Nframe);
            frame_noGI = RemoveGI(frame_parallel, Nfft*Nos, NGI*Nos);
            Frame_os = fft(frame_noGI,  Nfft*Nos);
            Frame = [Frame_os(Nfft/2+Nfft*(Nos-1)+1:Nfft*Nos, :); Frame_os(2:Nfft/2+1, :)];
            % Frame = fft(frame_parallel, Nfft);
%             for kkk = 1: Nframe
%                 Frame(:,kkk) = fftshift(Frame(:,kkk));
%             end
            Frame_demod = QPSKDemod(Frame, Nfft, Nframe);
            nobe = nobe + biterr(Frame_FDparallel(2:end,:),Frame_demod(2:end,:));

            % nobe = nobe + biterr(Frame_FDparallel(2:end,:),Frame_demod(2:end,:)); 
        end
        bers(kk, ii) = nobe/Maxiter/Nmod/Nk
        if ii == N_SNR
            for jjj = 1:N_Z
                CCDF(jjj) = sum(CF>Z_dBs(jjj))/Maxiter;
            end
            CCDF_CR(kk, :) = CCDF;
        end
    end
end
%% 画图
% figure(1)
% plot(abs(fft(frame_clipped)));
% figure(2)
% plot(abs(fft(frame_filter))); hold off
for kk = 1:N_CR
   gs = gss(kk);
   subplot(211), semilogy(Z_dBs,CCDF_CR(kk,:),[gs '-']), hold on
   subplot(212), semilogy(SNRdBs,bers(kk,:),[gs '-']), hold on
end 

2.代码分析

PAPR就是通频带上的概念，所以当然要涉及到通频带啦！我也就是为了这个写了这篇文章，可以说是为了这碗醋才包的这顿饺子

过采样是怎么完成的

我觉得本质上还是虚拟子载波，就像上面说的，原书中的仿真我也没有完全弄明白，我也是按照虚拟子载波的思路来处理的

各个输出的频域图像

因为我没有明白上载波的过程，所以我们画出来看看

如果我们无视所谓的fftshift呢

这里我取消了滤波，因为滤波了由于中心频率不对就已经错了

现在我们能明白fftshift的作用了吗！记得之前文章我们说过的吗，在fft中，所有的东西都是循环的，比如线性卷积变换为循环卷积。因此，想象一下我们正常的基带信号应该是正负半轴各一半，到了fft里负的一般自动循环位移到后面去啦。如果没有这一步，可以看到在通频带上信号的位置就错误了。

而如何从基带到通频带呢？这实际上就是把复包络变换到通带信号的变换

如何从通频带变换到基带呢？平移回来后滤波就好。

五. 剩下的一些疑问

1. 频域的横坐标应该如何画呢，我这样应该还是有问题的
2. dB量之间的关系？