OFDM发射机与接收机的matlab实现

clear;
close all;

EbN0 = [0:2:20];          %Eb/N0
N_iter = 1e5;
Nframe = 3;
Nused = 64;
Nbps = 4; M = 2^Nbps;     %调制阶数 = 2/4/6; QPSK/16-QAM/64-QAM
Nfft = 64;
Nsym = 64;
norms = [1 sqrt(2) 0 sqrt(10) 0 sqrt(42)];  %BPSK 4-QAM 16-QAM
Target_neb = 500;    %积累了一定的错误自动停止循环
sigPow = 0;
for i = 0:length(EbN0)
    randn('state',0);rand('state',0);
%     初始化错误比特数和总比特数
    Neb = 0;
    Ntb = 0;
%     迭代设计
    for m = 1:N_iter
        X = randi([0,M-1],1,Nused*Nframe);   %bit:整数向量
        Xmod = qammod(X,M,'gray')/norms(Nbps);
        
        kk1 = [1:Nused/2];          %前半有效载波序号
        kk2 = [Nused/2 + 1:Nused];  %后半有效载波序号
        kk3 = 1:Nfft;               %载波序号
        kk4 = 1:Nsym;               %符号周期序号
        
        for k = 1:Nframe
            X_shift  = [Xmod(kk2) Xmod(kk1)];
%         IFFT 调制
            x = ifft(X_shift);
            x_GI(kk4) = x;
            kk1 = kk1 + Nused;
            kk2 = kk2 + Nused;
            kk3 = kk3 + Nfft;
            kk4 = kk4 + Nsym;
        end
        
        %     信道
         y = x_GI;
        if i == 0   %只测量信号功率
            y1 = y(1:Nframe*Nsym);
            sigPow = sigPow + y1*y1';
            continue;
        end
%     y1 = x_GI(1:Nframe*Nsym);
   
%     sigPow = y1*y1';
    snr = EbN0(i) + 10*log10(Nbps*(Nused/Nfft));  %式(4.28)SNR vs.Eb/N0
%     noise_mag = 0;
    noise_mag = sqrt((10.^(-snr/10))*sigPow/2);
    
    y_GI = y + noise_mag * (randn(size(y))+1j*randn(size(y)));
%     Rx
    kk1 = 1:Nsym;
    kk2 = 1:Nfft;
    kk3 = 1:Nused;
    kk4 = Nused/2  + 1:Nfft;
    kk5 = [1:Nused/2];
    for k = 1:Nframe
%         FFT
    Y(kk2) = fft(y_GI(kk1));
    Y_shift = [Y(kk4) Y(kk5)];
    Xmod_r(kk3) = Y_shift;
    
    kk1 = kk1 + Nsym;
    kk2 = kk2 + Nfft;
    kk3 = kk3 + Nused;
    kk4 = kk4 + Nfft;
    kk5 = kk5 + Nfft;
    end
    X_r = qamdemod(Xmod_r*norms(Nbps),M,'gray');          %解调制
    Neb = Neb + sum(sum(de2bi(X_r,Nbps)~=de2bi(X,Nbps))); %计算误比特数
    Ntb = Ntb + Nused*Nframe * Nbps;                      %计算总比特数
    if Neb > Target_neb,break;end
    end
    
    if i == 0
        sigPow = sigPow/Nsym/Nframe/N_iter;
    else
        Ber = Neb/Ntb;
        fprintf('EbN0=%3d[dB],BER=%4d/%8d=%11.3e\n',EbN0(i),Neb,Ntb,Ber);
        Ber_buf(i) = Ber;
        if Ber < 1e-6 , break;end
    end
    
end
disp('Simulation finished')
EbN0dB = [0:1:20];
M = 2^Nbps;
ber_AWGN = ber_QAM(EbN0dB,M,'AWGN');
semilogy(EbN0dB,ber_AWGN,'r:'),hold on
semilogy(EbN0(1:i),Ber_buf,'b--s');
grid on;
axis([EbN0(1) EbN0(end) 1e-6 1])
legend('AWGN analytic','Simulation')
xlabel('EbN0[dB]'),ylabel('BER')
% plot(Ber_buf)