此程序是基于限幅滤波降低PAPR的程序。BER仿真一直出现问题。大佬们帮小弟看看好不

%CCDF_of_cliping_filtered_OFDM_signal.m
clear;
clf;
SNRdBs=(0:12);%SNR[dB]向量
N_SNR=length(SNRdBs);
Nblk=100;
CRs=(0.8:0.2:1.6);% 限幅比参数
N_CR=length(CRs);%限幅比个数
gss='*sdphv';%画图符号
b=2; %每个QAM符号的比特数
M=2^b; %和字符集大小
N=128; %FFT大小
Ncp=0; %CP大小（GI长度）
fs=1e6; %采样频率
L=8; %过采样因子
Tsym=(1/fs)*N; %OFDM符号周期
Ts=1/(fs*L); %OFDM采样周期
fc=2e6;
wc=2*pi*fc; %载波频率
t=(0:Ts:2*Tsym-Ts)/Tsym; %时间向量
A=modnorm(qammod((0:M-1),M),'avpow',1); %归一化因子  基于平均功率算出的比例因子
% mdmod = qammod(randi([0 M-1],1,N),M,'Gray');%QAM调制
% mddem = qamdemod(randi([0 M-1],1,N),M,'Gray');%QAM解调
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});%BPF系数

Clipped_errCnt=zeros(size(CRs));
ClippedFiltered_errCnt=zeros(size(CRs));
CF=zeros(1,Nblk);%没有限幅 滤波
CF_c=zeros(N_CR,Nblk);%只有限幅
CF_cf=zeros(N_CR,Nblk);%限幅加滤波
ber_analytic=berawgn(SNRdBs-10*log10(b),'qam',M);%返回QAM未编码的AWGN信道的BER
kk1=1:(N/2-Ncp)*L;%1-512
kk2=kk1(end)+1:N/2*L+N*L;%513-1536
kk3=kk2(end)+1:N*L/2+1536;%1537-2048
z=(2:0.1:16);
len_z=length(z);
%--------增加SNRdB进行循环--------%
for i=1:N_SNR%SNR[dB]向量(0-10)
%     msgbin=randi(b,1,N); %二进制序列%128*128

    SNRdB=SNRdBs(i);
    for ncf=0:2 %无限幅/限幅/限幅和滤波
        if ncf==2
            m=ceil(length(h)/2);%h是BPF的数值
        else  %ceil 是向离它最近的大整数圆整  返回大于或者等于指定表达式的最小整数
            m=1;
        end
        for cr=1:N_CR%限幅比个数
%             if ncf==0 & cr<1
%                 break;
%             end
            CR=CRs(cr);
            nobe=0;
            for nblk=1:Nblk %(i)
                msgbin=randi(b,1,N); %二进制序列%128*128
%            
                mdmod = qammod(msgbin,M,'Gray');%调制
%                 mddem = qamdemod(randi([0 M-1],1,N),M,'Gray');%解调
                X = A*modulate(mdmod,msgbin); %4QAM(QPSK)映射
                X(1)=0+1i*0; %不使用DC
                x=IFFT_oversampling(X,N,L);%过采样 IFFT
                x_b=add_CP(x,Ncp*L);%加CP =0
                x_b_os=[zeros(1,(N/2-Ncp)*L),x_b,zeros(1,N*L)];%前面和后面补零
                x_p=sqrt(2)*real(x_b_os.*exp(1i*wc*t));%数字上变换
                if ncf>0
                    x_p_c=clipping1(x_p,CR); %限幅
                    x_p=x_p_c;
                    if ncf>1
                        x_p_cf=ifft(fft(h,length(x_p)).*fft(x_p));%再加滤波
                        x_p=x_p_cf;
                    end
                end
                if i==N_SNR
                    CF(nblk)=PAPR(x_p);
                end
                %以下是用于解调
                y_p_n=[x_p(kk1) awgn(x_p(kk2),SNRdB,'measured') x_p(kk3)];%加噪声
                y_b=y_p_n.*exp(-1i*wc*t)/sqrt(2);%下变换
                Y_b=fft(y_b);
                y_b_z=ifft(zero_pasting(Y_b));
                y_b_t=y_b_z((N/2-Ncp)*L+m+[0:L:(N+Ncp)*L-1]);
                Y_b_f=fft(y_b_t(Ncp+1:end),N)*L;
                Y_b_bin=qamdemod(Y_b_f,M,'Gray');
                % Y_b_bin=demodulate(mddem,Y_b_f);%解调
                nobe=nobe+biterr(msgbin(:,2:end),Y_b_bin(:,2:end));%128*1 biterr函数只对无符号二进制数进行比较，即 0个数/1个数
%                 if ncf==0
%                     ber_no(i)=nobe/Nblk/(N-1)/b;
%                 elseif ncf==1
%                     ber_c(cr,i)=nobe/Nblk/(N-1)/b;
%                 else
%                     ber_cf(cr,i)=nobe/Nblk/(N-1)/b;
%                 end
            end

            if ncf==0
                ber_no(i)=nobe/Nblk/(N-1)/b;
            elseif ncf==1
                ber_c(cr,i)=nobe/Nblk/(N-1)/b;
            else
                ber_cf(cr,i)=nobe/Nblk/(N-1)/b;
            end
            if i==N_SNR
                for iz=1:len_z
                    CCDF(iz)=sum(CF>z(iz))/Nblk;
                end
                if ncf==0
                    CCDF_no=CCDF;
                    break;
                elseif ncf==1
                    CCDF_c(cr,:)=CCDF;
                else
                    CCDF_cf(cr,:)=CCDF;
                end
            end
        end
    end
end
figure(1)
semilogy(z,CCDF_no)
grid on,hold on
for cr=1:N_CR
    gs=gss(cr);
    semilogy(z,CCDF_c(cr,:),[gs '-'],z,CCDF_cf(cr,:),[gs ':'])%前实线  后虚线
    hold on
end

figure(2)
% semilogy(SNRdBs,ber_no,'y');
  semilogy(SNRdBs,ber_analytic,'r');
grid on,hold on
for cr=1:N_CR
    gs=gss(cr);
    semilogy(SNRdBs,ber_c(cr,:),[gs '-'],SNRdBs,ber_cf(cr,:),[gs ':'])
end