m基于PTS+TR的OFDM系统PAPR联合抑制算法matlab仿真

1.算法描述

       部分传输序列(Partial Transmit Sequence , PTS)由于其不受载波数量限制，并且能够有效的，无失真的降低OFDM信号峰均比，而受到广泛关注。部分传输序列算法(PTS)最初是由S.H.Muller和J.B.Huber于1997年提出。PTS算法的核心思想是将具有N个符号的输入序列按照一定的分割方式分割成V个子数据块，并且保持每个子数据块仍含有N个符号。然后对V个子数据块进行相位加权与合并处理，选择具有最小PAPR的一组符号进行传输，达到降低OFDM信号PAPR的目的。传统的PTS算法理论比较多，现成的资料也比较多，这里就不多做介绍了，通过仿真，对比PTS和没有PTS下。目前OFDM的PAPR主要算法有信号预畸变，信号扰码，编码三个方向来解决。

 

       在本课题中，我们将在传统PTS算法基础上引入了TR的思路到改进后的PTS算法中，引入的意义为：先预留出若干子载波来加载削峰信号，然后利用优化过的PTS算法对OFDM符号的PAPR进行抑制，之后再利用改进的TR算法对符号的PAPR进行进一步的抑制。整个算法的流程如下所示：

 

步骤一：加入门限，降低PTS算法的复杂度（但是这样会降低性能）

 

当满足要求：

 

算法就停止搜索，这样的话，就降低的算法的复杂度，但是会影响性能。

 

步骤二：加入限幅的方法

 

 

       通过这个方法，可以在步骤一的基础上，提高性能，使其在复杂度降低的前提下，保存系统的性能不变。

 

步骤三：改进PTS和TR的结合

 

    为了和TR结合，首先，PTS分组必须为随机分组，并随机的保留一定的预留子载波，然后先执行PTS，再执行TR。

 

步骤四：执行TR

 

       将得到的频域信号X进行IFFT变换得到时域信号x，对x的每个子载波上的数据限幅，对取反后的限幅差值进行N点FFT变换，得到的频域反向限幅差值信号的预留子载波上的数据即为削峰数据，用其替代X中预留子载波上的数据即可有效地消除峰值信号。

 

2.仿真效果预览

matlab2022a仿真结果如下：

3.MATLAB部分代码预览 `for k = 1:Nframes

    if mod(k,1000) == 0

       k/1000

    end  

    %产生数据源

    QPSK_Ind       = floor(length(Map_qpsk)*rand(1,Nfft))+1;

    %调制，这里为了研究PAPR性能，所以不加入编码模块和交织模块

    Qpsk_mod       = Map_qpsk(QPSK_Ind(1,:));            

    %随机分割

    tic;

    QPSK_Ind = randperm(Nfft);

    A        = zeros(1,Nfft);

    for v=1:Npts

        A(v,QPSK_Ind(v:Npts:Nfft)) = Qpsk_mod(QPSK_Ind(v:Npts:Nfft));

    end

    a       = ifft(A,[],2);

    %限幅

    [rr,cc] = size(a);

    for i = 1:rr

        for j = 1:cc

            if abs(a(i,j)) > Tho

               a(i,j) = Tho*(real(a(i,j)) + ij*imag(a(i,j)))/abs(a(i,j));

            end

        end

    end

    

    for n = 1:4^Npts

        %相位组合因子

        phase_temp        = Init_Phase(Data_back(n,:)).';

        if n == 1

           a_temp         = sum(a.*repmat(phase_temp,1,Nfft));

        else  

           a_temp         = a_temp + sum(a.*repmat(phase_temp,1,Nfft));

        end

        Signal_Power_temp = abs(a_temp.^2);

        Peak_Power_temp   = max(Signal_Power_temp,[],2);

        Mean_Power_temp   = mean(Signal_Power_temp,2);

        PAPR_temp         = 10*log10(Peak_Power_temp./Mean_Power_temp);

        if PAPR_temp < Th

           PAPR_pts(k) = PAPR_temp;

           X2          = a_temp;    

           break;

        end

    end

    %限幅

    [rr,cc] = size(X2);

    X2s     = X2;

    for i = 1:rr

        for j = 1:cc

            if abs(X2(i,j)) > Tho2

               X2s(i,j) = Tho2*(real(X2(i,j)) + ij*imag(X2(i,j)))/abs(X2(i,j));

            end

        end

    end

    X3 = X2s;

    

    Signal_Power_temp = abs(X3.^2);

    Peak_Power_temp   = max(Signal_Power_temp,[],2);

    Mean_Power_temp   = mean(Signal_Power_temp,2);

    PAPRs(k)          = 10*log10(Peak_Power_temp./Mean_Power_temp);

    times(k) = toc;

end

[cdf,PAPR] = ecdf(PAPRs);

figure;

semilogy(PAPR,1-cdf,'b','LineWidth',3);

xlabel('PAPR0[dB]');

ylabel('CCDF (Pr[PAPR>PAPR0])');

grid on;

title('有PAPR的时候的系统CCDF图');

save PAPR_Data_with_PAPR.mat PAPR cdf

%下面的代码是计算误码率的代码

Error    = zeros(1,length(SNR));

Rec      = zeros(1,Nfft);

PAPR_pts = zeros(1,min(Nframes,2000));

for ii = 1:length(SNR)

    Err_tmp = 0;

    for k=1:min(Nframes,2000)

%         RandStream.setDefaultStream(RandStream('mt19937ar','seed',k*ii));

        if mod(k,1000) == 0

           ii

           k/1000

        end

        %产生数据源

        QPSK_Dat     = floor(length(Map_qpsk)*rand(1,Nfft)) + 1;

        %调制，这里为了研究PAPR性能，所以不加入编码模块和交织模块

        Qpsk_mod     = Map_qpsk(QPSK_Dat);   

        

        %进行IFFT变换

        %随机分割

        QPSK_Ind = randperm(Nfft);

        A        = zeros(1,Nfft);

        for v=1:Npts

            A(v,QPSK_Ind(v:Npts:Nfft)) = Qpsk_mod(QPSK_Ind(v:Npts:Nfft));

        end

        a           = ifft(A,[],2);   

        %限幅

        [rr,cc] = size(a);

        for i = 1:rr

            for j = 1:cc

                if abs(a(i,j)) > Tho

                   a(i,j) = Tho*(real(a(i,j)) + ij*imag(a(i,j)))/abs(a(i,j));

                end

            end

        end

        

        for n = 1:4^Npts

            %相位组合因子

            phase_temp        = Init_Phase(Data_back(n,:)).';

            if n == 1

               a_temp         = sum(a.*repmat(phase_temp,1,Nfft));

            else  

               a_temp         = a_temp + sum(a.*repmat(phase_temp,1,Nfft));

            end

            Signal_Power_temp = abs(a_temp.^2);

            Peak_Power_temp   = max(Signal_Power_temp,[],2);

            Mean_Power_temp   = mean(Signal_Power_temp,2);

            PAPR_temp         = 10*log10(Peak_Power_temp./Mean_Power_temp);

            if PAPR_temp < Th

               PAPR_pts(k)    = PAPR_temp;

               X2             = a_temp;    

               break;

            end

        end

        %限幅

        [rr,cc] = size(X2);

        X2s     = X2;

        for i = 1:rr

            for j = 1:cc

                if abs(X2(i,j)) > Tho2

                   X2s(i,j) = Tho2*(real(X2(i,j)) + ij*imag(X2(i,j)))/abs(X2(i,j));

                end

            end

        end

        X3 = X2s;

 

        R = X3;

        %通过高斯信道

        Dat_Ifft     = awgn(R,SNR(ii),'measured');

        %模拟实际的接收端的畸变

        Dat_Ifft2    = Dat_Ifft;

        if PAPR_pts(k) > 8+Tho+Tho2%瞬时功率过大，则畸变

           Dat_Ifft2 = randn(1,Nfft) + ij*randn(1,Nfft);

        end

        

        %fft变换

        Dat_fft      = fft(Dat_Ifft2,[],2);

        %解调

        I            = sign(real(Dat_fft)).*(abs(real(Dat_fft))>0.5);

        Q            = sign(imag(Dat_fft)).*(abs(imag(Dat_fft))>0.5);

        for i = 1:Nfft

            if I(i) ==  1 & Q(i) ==  0

               Rec(i) = 1;

            end

            if I(i) == -1 & Q(i) ==  0

               Rec(i) = 2;

            end

            if I(i) == 0  & Q(i) ==  1

               Rec(i) = 3;

            end

            if I(i) == 0  & Q(i) == -1

               Rec(i) = 4;

            end            

        end

    Err_tmp = Err_tmp + length(find(QPSK_Dat~=Rec));    

    end

    Error(ii) = Err_tmp/min(Nframes,2000)/Nfft;

end

01_060_m`