ofdm_awgn.m

文件夹“AWGN”说明了系统误码率在均衡前后随SNR增加的变化。 （1）绘制误码率曲线。运行m文件：main_picture,调用主函数：ofdm_AWGN。 （2）绘制星座点图。运行m文件：sn

%      !------------------Fixed gain-------!
%      !                                   !
%s(t)--!signal in                         Sum--------r(t)signal out
%      !                                   !
%      !---------Delay1----Variable gain----!
%      !                   or attenuation
%      !---------Delay2----Variable gain----!
%                         or attenuation
 %clear;clc
 %nl=3;
 %erros_curve=zeros(3,6400);

% SNRDB=[-2:2]
%erro_un=zeros(1,length(SNRDB));
% erro_eq=zeros(1,length(SNRDB));

function [BER,BER_EQU] = ofdm_AWGN(SNR_DB)
lll=1;
for SNR_dB = 0:20;
    
%********************** preparation part ***************************
bitrate = 24e6;         % 系统传输速率
TC = 3.2e-6;            % Duration of one OFDM symbol,1/TC is the carrier frequency
TP = 800e-9;            % Duration of the cyclic prefix [S]
TG = 800e-9;            % Tolal guard time [S]
nf = 10;               % 一帧100个符号
ml =4;                  % 映射类别:BPSK=1,QPSK=2,16QAM=4,64QAM=6 
over=4;                 % 过采样倍数
NFFT = 128;             % FFT点数
A = 1;                  % Amplitude of the rectangular impulse response [V]
fp = 1e8;               % 载频
fs =40e8;              % 采样频率
ts = 1/fs;              % 采样周期
snrdb= SNR_dB;                % 高斯白噪声的信噪比
noise_choice=0;       % 乘性干扰的种类：1——复乘性干扰，0——实乘性干扰
numTaps =  32 ;		% channel order
Mu = 0.03;          % iteration step size
%--------------------------------------------------------------------------
eob = 0;    % Number of error bit
eob_equ=0;
nob = 0;    % Number of transmitted bit
eop=0;      % Number of error packet
nop=0;      % Number of transmitted packet
%--------------------------------------------------------------------------
for loop1=1:10
%--------------------------------------------------------------------------
nb = bitrate*(TC+TP);   % 一个符号的bit数
numcarrier = nb/ml;     % 子载波数
numbits = nb*nf;        % 一帧的bit数
%--------------------------------------------------------------------------
tc = TC / NFFT;            % Duration per chip
ncp = floor(TP/tc);        % Number of tones of the cyclic prefix
n = (-ncp+1:1:NFFT);     
N = length(n);             % Total number of tones per symbol
%--------------------------------------------------------------------------
if ncp>0
   tc = (TC+TP)/N;         % Tone duration
end
chipsample = floor(tc*fs);      % 每个切普的采样点数
chipresid = floor((TG-TP)*fs);  % Samples of the residual part
symbolsamp = (chipsample*N)+chipresid; % 一个符号的采样点数
framesamp = symbolsamp * nf;           % 一帧的采样点总数
%************************** transmitter ***********************************

%**************************************************************************
%                            信号源
%**************************************************************************
%[ge_bits] = rand_bits(numbits);
ge_bits=randint(1,numbits); % 产生二进制数据流
% ge_bits=zeros(1,numbits);
%ge_bits=ones(1,numbits);
%g=[0 0 0 0 0 1 1 0 1 0 0 1 ] ;
%ge_bits=repmat(g,1,numbits/length(g));
%--------------------------------------------------------------------------
%                          画图部分 --信号源
%--------------------------------------------------------------------------
%pause 
%figure
%t1 = linspace (0,nf*(TC+TP),length(ge_bits));
%stairs(t1,ge_bits)                 %图1 
% ylim([-0.3 1.5])
% X = xlabel('Time [S]');
% Y = ylabel('Amplitude [V]');
% title('发送的信息')
%****************************************************************
%                          映射
%****************************************************************
S =yinshe(ge_bits,ml);
%--------------------------------------------------------------------------
%                             画图部分 --映射
%--------------------------------------------------------------------------
%pause
%figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 星座图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot(real(S),imag(S),'.')
% xlim([-1.5 1.5])
% ylim([-1.5 1.5])

%t2 = linspace (0,nf*(TC+TP),length(S));
%stairs(t2,abs(S));
%ylim([0 2])
%X = xlabel('Time [S]');
%Y = ylabel('Amplitude [V]');
%title( '映射后幅度')

%pause
%figure
%plot(t2,angle(S))
%ylim([-3 4.2])
%X = xlabel('Time [S]');
%Y = ylabel('Phase[rad]');
%title( '映射后相位')
%**************************************************************************
%                            串并转换
%**************************************************************************
para_S = reshape(S,numcarrier,nf);  % Serial to Parallel Conversion
%**************************************************************************
%                           插0使长度变成128
%**************************************************************************
para_S0 = zeros(NFFT,nf);
up1 = floor(numcarrier/2);
down1 = numcarrier-up1;
para_S0(1:up1,:) = para_S(1:up1,:);
para_S0(NFFT-down1+1:NFFT,:) = para_S(numcarrier-down1+1:numcarrier,:); 
%**************************************************************************
%                              4倍过采样
%**************************************************************************
 up2 = floor(NFFT/2);
 down2 = NFFT-up2;
 FS = over*NFFT;                          %FS=512
 FS_cp =over*(ncp+NFFT);                  %FS_cp=640
 ifftlen = nf*FS_cp;                     %ifftlen=6400
 ifft_data = [zeros(1,ifftlen)'];
 for b = 1:nf                             %循环说明一个一个发送符号
    ch = para_S0(:,b); 
    Czp = zeros(FS,1);                  %4倍过采样中间补0
    Czp(1:up2)=ch(1:up2);
    Czp(FS-down2+1:FS)=ch(NFFT-down2+1:NFFT);
 %*************************************************************************
 %                                 IFFT
 %*************************************************************************
    C = ifft(Czp,FS);     % IFFT of the zero-padded input
    C = sqrt(FS)*C;
%**************************************************************************
%                            插入CP
%**************************************************************************
     if ncp>0 
        C_cp=zeros(length(C)+over*ncp,1);
        C_cp(1:(over*ncp))=C(over*NFFT-(over*ncp)+1:over*NFFT);  
        C_cp(over*ncp+1:length(C_cp))=C;
     else
        C_cp=C;
     end
    ifft_data((1+(b-1)*FS_cp):(FS_cp+(b-1)*FS_cp)) = C_cp(:);
% end

%**************************************************************************
%                           并串转换
%**************************************************************************
XM=ifft_data((1+(b-1)*FS_cp):(FS_cp+(b-1)*FS_cp))';
%--------------------------------------------------------------------------
%                        画图程序--- IFFT后
%--------------------------------------------------------------------------
% t3 = linspace (0,ts*ifftlen,length(XM));
% subplot(2,1,1)
% plot(t3,abs(XM))                                      %图3
% X = xlabel('Time [S]');
% Y = ylabel('Amplitude [V]');
% title('IFFT 后模值')
% 
% subplot(2,1,2)
% plot(t3,angle(XM))                                      %图3
% X = xlabel('Time [S]');
% Y = ylabel('Phase [V]');
% title('IFFT 后相位')
%************************************************************************** 
%                             上变频
%**************************************************************************
[Stx,time] = modulation(XM,ts,ifftlen,fp);

%**************************************************************************
%                             四径信道 
%**************************************************************************
 signal_in=Stx;
 %signal_out=multpath(Stx,add_delay2,add_delay3,add_delay4,noise_noise_choice);
 for ii = 1:length(signal_in) 
   g2(ii)=0.6*rand(1);
   g3(ii)=0.6*rand(1);
   g4(ii)=0.6*rand(1);
   g5(ii)=0.6*rand(1);
   g6(ii)=0.6*rand(1);
   g7(ii)=0.6*rand(1);
   signal2(ii)=(g2(ii)+i*g3(ii)*noise_choice)*signal_in(ii);
   signal3(ii)=(g4(ii)+i*g5(ii)*noise_choice)*signal_in(ii);
   %signal4(ii)=(g6(ii)+i*g7(ii)*noise_choice)*signal_in(ii);
end    
 tau2=round(10*abs(randn(1))+154);        % variable delay(phase shift)
  %tau3=round(10*abs(randn(1))+161);
 tau3=round(10*abs(randn(1))+172);
% tau2=128+loop1;
% tau3=138+2*loop1;
 %signal_out=signal_in+[signal2(end-tau2+1:end) signal2(1:end-tau2)]+[signal3(end-tau3+1:end) signal3(1:end-tau3)];
 signal_out=signal_in;
   %signal_out=signal_in+signal2+signal3;
%**************************************************************************
%                             加上高斯白噪声
%**************************************************************************
 Ec = 1;
 ecn = 10.^(snrdb/10);
 noise_var = Ec/ecn;
 noise_var = sqrt(noise_var)/2;
%[Rx_data] = add_whitenoise(signal_out,noise_var);
 Rx_data = awgn(signal_out,snrdb);
signal_out(:)=Rx_data(:);
%**************************************************************************
%                             Receiver  
%**************************************************************************

%********************************************************************
%                         下变频
%********************************************************************
%[R_XM] = demod(Stx,ts,ifftlen,fp);
[R_XM] = demod(signal_out,ts,ifftlen,fp);
%[iR_XM,qR_XM] = demod(Stx,ts,ifftlen,fp);
%********************************************************************
%                         低通滤波器
%********************************************************************
%[iRx_XM,qRx_XM] =lowpass(iR_XM,qR_XM,fp,fs);
Rx_XM((1+(b-1)*FS_cp):(FS_cp+(b-1)*FS_cp))=R_XM(:);
 end 
%**************************************************************************
%                         把nf个符号按顺序排列
%**************************************************************************
Rx_X = Rx_XM';                           %理论上,Rx_X=ifft_data;
%**************************************************************************
%**************************************************************************
%                           串并转换
%**************************************************************************
nc = floor(length(Rx_X)/FS_cp);
Rx_X = reshape(Rx_X,FS_cp,nc);
Rx_ch = zeros(NFFT,nc);
%*************************************************************************** 
%                          去循环前缀 
%**************************************************************************
for m = 1:nc
   Rx_C_cp = Rx_X(:,m);
   if ncp>0
        Rx_C = zeros(over*NFFT,1);
        Rx_C(:) = Rx_C_cp(over*ncp+1:length(Rx_C_cp));
    else
        Rx_C = Rx_C_cp;
 end
%********************************************************************   
%                             FFT 
%********************************************************************
    %FS = 4*NFFT;
    Rx_C = Rx_C/sqrt(FS);
    Rx_Czp = fft(Rx_C,FS);
    %disp(Rx_Czp);
%**************************************************************************
%                              还原过采样
%**************************************************************************
  % up2 = floor(NFFT/2);
  % down2 = NFFT- up2;
    Rx_c0 = zeros(NFFT,1);
    Rx_c0(1:up2) = Rx_Czp(1:up2);
    Rx_c0(up2+1:NFFT) = Rx_Czp(FS-down2+1:FS); 
    Rx_ch(:,m)=Rx_c0(:);
end
%**************************************************************************
%                         去掉128中间的0
%**************************************************************************
    Rx_S0 = zeros(numcarrier,nf);
    %up1 = floor(numcarrier/2);
    %down1 = numcarrier-up1;
    Rx_S0(1:up1,:) = Rx_ch(1:up1,:);
    Rx_S0(up1+1:numcarrier,:) = Rx_ch(NFFT-down1+1:NFFT,:);
%**************************************************************************
%                               并串转换
%**************************************************************************
Rx_S = reshape(Rx_S0,1,numcarrier*nc);
%Rx_S = sign(Rx_S);
%--------------------------------------------------------------------------
%                         FFT后画图程序