echo off;clear all;close all;clc;
fprintf( 'OFDM仿真\n') ;
tic
% --------------------------------------------- %
% 参数定义 %
% --------------------------------------------- %
% Initialize the parameters
NumLoop = 100;
NumSubc =600;
SyncDelay = 0;
N=1024;
T=0.2; %码元宽度
T_guard=0.05; %保护间隔
NumCP = 2*N*T_guard/T;
%-----------------------------------------------
%频带范围:2K-8K;载波数为600;频率间隔10HZ;
% 子载波数 600
% 位数/ 符号 2
% 符号数/ 载波 100
% 训练符号数 0
%做ifft/fft的数据为N=1024,最接近600的2的N次倍
% 循环前缀长度 NumCP = 2*N*T_guard/T;
% 调制方式 4-QAM
% 多径信道数 5
%IFFT Size 1024
%-----------------------------------------------
% --------------------------------------------- %
% QAM MODULATION %
% --------------------------------------------- %
% Generate the random binary stream for transmit test
BitsTx = floor(rand(1,NumLoop*NumSubc*2)*2);%产生0,1随机数
% Modulate (Generates QAM symbols)
% input: BitsTx(1,NumLoop*NumSubc); output:SymQAM(NumLoop,NumSubc/2)
SymQAMtmp = reshape(BitsTx,2,NumLoop*NumSubc).';
SymQAMtmptmp = bi2de(SymQAMtmp,2,'left-msb');
%--------------------------------------------------------------------------
% 函数说明:
% bin2dec(binarystr) interprets the binary string binarystr and returns the equivalent decimal number.
% bi2de是把列向量的每一个元素都由2进制变为10进制
% D = BI2DE(...,MSBFLAG) uses MSBFLAG to determine the input orientation.
% MSBFLAG has two possible values, 'right-msb' and 'left-msb'. Giving a 'right-msb' MSBFLAG does not change the function's default behavior.
% Giving a 'left-msb' MSBFLAG flips the input orientation such that the
% MSB is on the left.
% D = BI2DE(...,P) converts a base P vector to a decimal value.
% Examples:
%>> B = [0 0 1 1; 1 0 1 0];
%>> T = [0 1 1; 2 1 0];
% >> D = bi2de(B) >> D = bi2de(B,'left-msb') >> D = bi2de(T,3)
% D = D = D =
% % 12 3 12
% % 5 10 5
%--------------------------------------------------------------------------
% QAM modulation
% 00->-1-i,01->-1+i,10->1-i,11->1+i
% 利用查表法进行QAM星座映射
QAMTable = [-1-i -1+i 1-i 1+i];
SymQAM = QAMTable(SymQAMtmptmp+1);
% --------------------------------------------- %
% IFFT
% --------------------------------------------- %
% input: SymQAM(NumLoop,NumSubc); output: SymIFFT(NumSubc,NumLoop)
SymIFFT = zeros(1024*2,NumLoop);%?
SymIFFTtmp = reshape(SymQAM,NumSubc,NumLoop);
T_s=zeros(1,NumSubc);
T_s=SymIFFTtmp(:,1);%添加训练序列
SymIFFTtmptmp = zeros(1024*2,NumLoop);%?
SymIFFTtmptmp(1:200,:) = SymIFFT(1:200,:); % 实数
SymIFFTtmptmp(201:800,:) = SymIFFTtmp; % 实数
SymIFFTtmptmp(801:1024,:)=SymIFFT(801:1024,:);
% 构造了一个2*N by NumLoop点的矩阵,这么安排矩阵的目的是为了构造共轭对称矩阵;
% 共轭对称矩阵的特点是 在ifft/fft的矢量上2*N点的矢量
% 在1,N+1点必须是实数
% 2至N点 与 N+2至2*N点关于N+1共轭对称
SymIFFTtmptmp(N+1,:) = SymIFFT(N+1,:);
SymIFFTtmptmp(N+2:N*2,:) = flipdim(conj(SymIFFTtmptmp(2:N,:)),1);
%--------------------------------------------------------------------------
% 函数说明:
% B = flipdim(A,dim) returns A with dimension dim flipped.
% When the value of dim is 1, the array is flipped row-wise down. When dim is 2,
% the array is flipped columnwise left to right. flipdim(A,1) is the same as
% flipud(A), and flipdim(A,2) is the same as fliplr(A).
%--------------------------------------------------------------------------
% >> a = [1 2 3; 4 5 6; 7 8 9; 10 11 12]
% a =
% 1 2 3
% 4 5 6
% 7 8 9
% 10 11 12
% >> b = flipdim(a,1)
% b =
% 10 11 12
% 7 8 9
% 4 5 6
% 1 2 3
SymIFFT = ifft(SymIFFTtmptmp,N*2,1);
% --------------------------------------------- %
% Add cyclic prefix %
% --------------------------------------------- %
% input: SymIFFT(2*N,NumLoop); output: SymCP(2*N + NumCP,NumLoop)
%plot(SymIFFT)
NumAddPrefix = 2*N + NumCP;
SymCP = zeros(NumAddPrefix,NumLoop);
RowPrefix = (2*N - NumCP + 1):2*N;
SymCP = [SymIFFT(RowPrefix,:);SymIFFT];%将时域波形加入保护间隔
% --------------------------------------------- %
% Go through the channel %
% --------------------------------------------- %
% input: SymCP(2*N + NumCP,NumLoop); output: SymCh(1,(2*N + NumCP)*NumLoop)
%SymChtmp = zeros(1,(2*N + NumCP)*NumLoop);
SymChtmp = SymCP(:).'; % 进行这个转置操作之后就成了一个矢量
% 相当于把矩阵的列向量依次排列 改变为一个行向量,100行,2*1024列
a1=0.83574;a2=1;a3=0.48634;a4=0.27974;a5=0.26422;
t1=0;t2=0.00983;t3=0.01640;t4=0.02608;t5=0.03117;
m1=round(t1*2*N/T);m2=round(t2*2*N/T);
m3=round(t3*2*N/T);
m4=round(t4*N*2/T);
m5=round(t5*N*2/T);
Ch = [a1,zeros(1,m2-1),a2,zeros(1,m3-m2-1),a3,zeros(1,m4-m3-1),a4,zeros(1,m5-m4-1),a5];
%Ch=[1];
SymChtmptmp = filter(Ch,1,SymChtmp);
%--------------------------------------------------------------------------
% 函数说明:
% Firlter data with an infinite impulse response (IIR) or finite impulse response
% (FIR) filter
% y = filter(b,a,X) filters the data in vector X with the filter described by
% numerator coefficient vector b and denominator coefficient vector a. If a(1) is
% not equal to 1, filter normalizes the filter coefficients by a(1). If a(1) equals
% 0, filter returns an error.
%--------------------------------------------------------------------------
% If X is a matrix, filter operates on the columns of X. If X is a multidimensional
% array, filter operates on the first nonsingleton dimension.
%-------------------------------------------------------------------------
% Add the AWGN
BerSnrTable = zeros(20,3);
s1=zeros(1,20);
s2=zeros(1,20);
for snr=0:19; % = SNR + 10*log10(log2(2));
BerSnrTable(snr+1,1) = snr;
SymCh = awgn(SymChtmptmp,snr,'measured');
%--------------------------------------------------------------------------
% 函数说明:
% AWGN Add white Gaussian noise to a signal.
% Y = AWGN(X,SNR) adds white Gaussian noise to X. The SNR is in dB.
% The power of X is assumed to be 0 dBW. If X is complex, then
% AWGN adds complex noise.
% ------------------------------------------------------------------------
% Y = AWGN(X,SNR,SIGPOWER) when SIGPOWER is numeric, it represents
% the signal power in dBW. When SIGPOWER is 'measured', AWGN measures
% the signal power before adding noise.
% -------------------------------------------------------------------------
% Y = AWGN(X,SNR,SIGPOWER,STATE) resets the state of RANDN to STATE.
%
% Y = AWGN(..., POWERTYPE) specifies the units of SNR and SIGPOWER.
% POWERTYPE can be 'db' or 'linear'. If POWERTYPE is 'db', then SNR
% is measured in dB and SIGPOWER is measured in dBW. If POWERTYPE is
% 'linear', then SNR is measured as a ratio and SIGPOWER is measured
% in Watts.
%
% Example: To specify the power of X to be 0 dBW and add noise to produce an SNR of 10dB, use:
% X = sqrt(2)*sin(0:pi/8:6*pi);
% Y = AWGN(X,10,0);
% Example: To specify the power of X to be 0 dBW, set RANDN to the 1234th
% state and add noise to produce an SNR of 10dB, use:
% X = sqrt(2)*sin(0:pi/8:6*pi);
% Y = AWGN(X,10,0,1234);
% Example: To specify the power of X to be 3 Watts and add noise to
% produce a linear SNR of 4, use:
% X = sqrt(2)*sin(0:pi/8:6*pi);
% Y = AWGN(X,4,3,'linear');
% Example: To cause AWGN to measure the power of X, set RANDN to the
% 1234th state and add noise to