www.pudn.com > MIMO-OFDM(simulinkANDmatlab).rar > ofdm.m, change:2007-12-24,size:14591b


echo off;clear all;close all;clc; 
fprintf( 'OFDM仿真\n') ; 
tic 
% --------------------------------------------- % 
%                   参数定义                     % 
% --------------------------------------------- % 
% Initialize the parameters 
NumLoop = 1000; 
NumSubc = 128; 
NumCP = 8; 
SyncDelay = 0; 
% 子载波数            128 
% 位数/ 符号          2 
% 符号数/ 载波        1000 
% 训练符号数          0 
% 循环前缀长度        8   (1/16)*T 
% 调制方式            4-QAM 
% 多径信道数          3  
% IFFT Size           128  
% 信道最大时延        2 
% --------------------------------------------- % 
%                   QAM MODULATION              % 
% --------------------------------------------- % 
% Generate the random binary stream for transmit test 
BitsTx = floor(rand(1,NumLoop*NumSubc)*2); 
% Modulate (Generates QAM symbols) 
% input: BitsTx(1,NumLoop*NumSubc); output: SymQAM(NumLoop,NumSubc/2) 
SymQAMtmp = reshape(BitsTx,2,NumLoop*NumSubc/2).'; 
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/2); output: SymIFFT(NumSubc,NumLoop) 
SymIFFT = zeros(NumSubc,NumLoop); 
SymIFFTtmp = reshape(SymQAM,NumSubc/2,NumLoop); 
SymIFFTtmptmp = zeros(NumSubc,NumLoop); 
SymIFFTtmptmp(1,:) = real(SymIFFTtmp(1,:)); % 实数 
SymIFFTtmptmp(NumSubc/2+1,:) = imag(SymIFFTtmp(1,:)); % 实数 
% 这么安排矩阵的目的是为了构造共轭对称矩阵 
% 共轭对称矩阵的特点是 在ifft/fft的矢量上 N点的矢量 
% 在0,N/2点必须是实数 一般选为0 
% 1至N/2点 与 (N/2)+1至N-1点关于N/2共轭对称 
SymIFFTtmptmp(2:NumSubc/2,:) = SymIFFTtmp(2:NumSubc/2,:); 
SymIFFTtmptmp((NumSubc/2+2):NumSubc,:) = flipdim(conj(SymIFFTtmp(2:NumSubc/2,:)),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,NumSubc,1); 
% --------------------------------------------- % 
%             Add cyclic prefix                 % 
% --------------------------------------------- % 
% input: SymIFFT(NumSubc,NumLoop); output: SymCP(NumSubc + NumCP,NumLoop) 
NumAddPrefix = NumSubc + NumCP; 
SymCP = zeros(NumAddPrefix,NumLoop); 
RowPrefix = (NumSubc - NumCP + 1):NumSubc; 
SymCP = [SymIFFT(RowPrefix,:);SymIFFT]; 
% --------------------------------------------- % 
%             Go through the channel            % 
% --------------------------------------------- % 
% input: SymCP(NumSubc + NumCP,NumLoop); output: SymCh(1,(NumSubc + NumCP)*NumLoop) 
SymCh = zeros(1,(NumSubc + NumCP)*NumLoop); 
SymChtmp = SymCP(:).';   % 进行这个转置操作之后就成了一个矢量 
                      % 相当于把矩阵的列向量依次排列 改变为一个行向量 
Ch = [1 1/2 1/4]; 
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); 
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 produce a linear SNR of 4, use: 
%              X = sqrt(2)*sin(0:pi/8:6*pi); 
%              Y = AWGN(X,4,'measured',1234,'linear'); 
% --------------------------------------------- % 
%            Remove Guard Intervals             % 
% --------------------------------------------- % 
% input: SymCh(1,(NumSubc + NumCP)*NumLoop); output: SymDeCP(NumSubc,NumLoop) 
SymDeCP = zeros(NumSubc,NumLoop); 
SymDeCPtmp = reshape(SymCh,NumSubc + NumCP,NumLoop); 
SymDeCP = SymDeCPtmp((NumCP+1+SyncDelay):NumAddPrefix+SyncDelay,:); 
% --------------------------------------------- % 
%                     FFT                       % 
% --------------------------------------------- % 
% input: SymDeCP(NumSubc,NumLoop); output: SymFFT(NumSubc,NumLoop) 
SymFFT = fft(SymDeCP,NumSubc,1); 
% --------------------------------------------- % 
%        Make Decision(Include DeQAM)           % 
% --------------------------------------------- % 
% SymFFT(NumSubc,NumLoop); output: SymDec(NumSubc,NumLoop) 
SymDec = zeros(NumSubc,NumLoop); 
SymEqtmp(1,:) = SymFFT(1,:)+i*SymFFT(NumSubc/2+1,:); 
SymEqtmp(2:NumSubc/2,:) = SymFFT(2:NumSubc/2,:); 
for m = 1:NumLoop 
    for n = 1:NumSubc/2 
        Real = real(SymEqtmp(n,m)); 
        Imag = imag(SymEqtmp(n,m)); 
         
        if( abs((Real -1)) < abs((Real +1))) 
            SymDec(2*n-1,m) = 1; 
        else 
            SymDec(2*n-1,m) = 0; 
        end 
        if( abs((Imag -1)) < abs((Imag +1  )) )   
            SymDec(2*n,m) = 1; 
        else 
            SymDec(2*n,m) = 0; 
        end 
    end 
end 
% --------------------------------------------------------------------- 
% Test by lavabin 
% Another way to DeQAM 
%  QAMTable = [-1-i -1+i 1-i 1+i]; 
%  00->-1-i,01->-1+i,10->1-i,11->1+i 
TestSymDec = zeros(NumSubc,NumLoop); 
TestSymEqtmp(1,:) = SymFFT(1,:)+i*SymFFT(NumSubc/2+1,:); 
TestSymEqtmp(2:NumSubc/2,:) = SymFFT(2:NumSubc/2,:); 
TestSymEqtmp1 = reshape(TestSymEqtmp,1,NumSubc*NumLoop/2); 
min_d = zeros(size(TestSymEqtmp1)); 
min_ddd = zeros(1,NumSubc*NumLoop); 
d = zeros(4,1); 
min_index = 0; 
for ii = 1:1:(NumSubc*NumLoop/2) 
    for jj = 1:4 
        d(jj) = abs(TestSymEqtmp(ii) - QAMTable(jj)); 
    end 
     [min_d(ii),min_index] = min(d); 
% % [Y,I] = MIN(X) returns the indices of the minimum values in vector I.     
      switch min_index 
  case 1 
     min_ddd(2*ii-1) = 0 ; 
         min_ddd(2*ii) =  0 ; 
  case 2 
     min_ddd(2*ii-1) = 0 ; 
         min_ddd(2*ii) =  1 ; 
  case 3 
     min_ddd(2*ii-1) = 1 ; 
         min_ddd(2*ii) =  0 ; 
  case 4 
     min_ddd(2*ii-1) = 1 ; 
         min_ddd(2*ii) =  1 ; 
          otherwise 
      fprintf('Impossible error!!! \n\n'); 
      end 
end 
%-------------------------------------------------------------------- 
% 函数说明: 
% % C = min(A) returns the smallest elements along different dimensions of an 
% % array. 
% % If A is a vector, min(A) returns the smallest element in A. 
% % If A is a matrix, min(A) treats the columns of A as vectors, returning a row 
% % vector containing the minimum element from each column. 
% % [C,I] = min(...) finds the indices of the minimum values of A, and returns 
% % them in output vector I. If there are several identical minimum values, the 
% % index of the first one found is returned. 
% Bit Error 
BitsRx = zeros(1,NumSubc*NumLoop); 
BitsRx = SymDec(:).'; 
[Num,Ber] = symerr(BitsTx,BitsRx) 
BerSnrTable(snr+1,2) = Num ; 
BerSnrTable(snr+1,3) = Ber ; 
end 
%-------------------------------------------------------------------- 
% Test by lavabin 
if min_ddd == BitsRx  
    fprintf('DeQAM two ways the same results \n\n'); 
else 
     fprintf('DeQAM two ways the different results'); 
end 
%-------------------------------------------------------------------- 
figure(1); 
subplot(2,1,1); 
semilogy(BerSnrTable(:,1),BerSnrTable(:,2),'o-'); 
subplot(2,1,2); 
semilogy(BerSnrTable(:,1),BerSnrTable(:,3),'o-'); 
%-------------------------------------------------------------------- 
time_of_sim = toc 
echo on; 
% --------------------------------------------- % 
%                   The END                     % 
% --------------------------------------------- % 
%    本程序中得到的收端OFDM信号的频谱波形,是与其发端信号的排步有关的。在发端的 
% 载波安排上,128个载波有前后各32个载波是null载波(如果这前后各32各载波是带外频段, 
% 那么理论上它们都应该是零!),中间的64个载波是数据载波。这样的排步很明显就是一个 
% 两边低,中间高的频谱形式。所以,收端也应该是这个轮廓。 
clc;clear all;close all;echo off;tic; 
% ------------------------------------------------------------------- 
%                      Parameter Definition 
% -------------------------------------------------------------- 
Fd    = 1;             % symbol rate (1Hz) 
Fs    = 1*Fd;          % number of sample per symbol 
M     = 4;             % kind(range) of symbol (0,1,2,3) 
Ndata = 1024;          % all transmitted data symbol  
Sdata = 64;            % 64 data symbol per frame to ifft 
Slen  = 128;           % 128 length symbol for IFFT  
Nsym  = Ndata/Sdata;   % number of frames -> Nsym frame 
GIlen = 144;           % symbol with GI insertion  GIlen = Slen + GI 
GI    = 16;            % guard interval length 
% ---------------------------------------------------------------- 
%                        Vector Initialization 
% ---------------------------------------------------------------- 
X  = zeros(Ndata,1); 
Y1 = zeros(Ndata,1); 
Y2 = zeros(Ndata,1); 
Y3 = zeros(Slen,1); 
z0 = zeros(Slen,1); 
z1 = zeros(Ndata/Sdata*Slen,1); 
g  = zeros(GIlen,1); 
z2 = zeros(GIlen*Nsym,1); 
z3 = zeros(GIlen*Nsym,1); 
% random integer generation by M kinds  
X = randint(Ndata, 1, M); 
% digital symbol mapped as analog symbol 
% Y1 is a Ndata-by-2 matrix, is changed into Y2 by "amodce" 
Y1 = modmap(X, Fd, Fs, 'qask', M); 
% covert to complex number 
Y2 = amodce(Y1,1,'qam'); 
% figure(1); 
% scatterplot(Y2,length(Y2),0,'bo');grid on; 
scatterplot(Y2,Fd,0,'bo');grid on; 
title('4-QAM Constellation'); 
Tx_spectrum = zeros(size(Y3)); 
for j=1:Nsym; 
    for i=1:Sdata; 
        Y3(i+Slen/2-Sdata/2,1)=Y2(i+(j-1)*Sdata,1); 
        Tx_spectrum = Tx_spectrum + abs(Y3); 
    end 
    z0=ifft(Y3); 
     
    for i=1:Slen;    % generate time-domain vector, z1, without GI 
        z1(((j-1)*Slen)+i)=z0(i,1); 
    end 
    % 
    for i=1:Slen; 
        g(i+16)=z0(i,1); 
    end 
     
    for i=1:GI; 
        g(i)=z0(i+Slen-GI,1);     
    end 
    for i=1:GIlen;  % generate time-domain vector, z2, with GI 
        z2(((j-1)*GIlen)+i)=g(i,1); 
    end 
     
end 
% graph on time domain 
figure(2); 
f = linspace(-Sdata,Sdata,length(z1)); 
plot(f,abs(z1)); 
    
Y4 = fft(z1);  % FFT operation, at receiver 
% if Y4 is under 0.01 Y4=0.001 
for j=1:Ndata/Sdata*Slen; 
if abs(Y4(j)) < 0.01 
    Y4(j)=0.01; 
end 
end 
Y4 = 10*log10(abs(Y4));  %  Y4 is used for spectrum display. 
% graph on frequency domain 
figure(3); 
f = linspace(-Sdata,Sdata,length(Y4)); 
plot(f,Y4);grid on; 
axis([-Slen/2 Slen/2 -20 20]); 
title('Received OFDM signal spectrum'); 
figure(4); 
f = linspace(-Sdata,Sdata,length(Y3)); 
plot(f,abs(Y3),'b-','LineWidth',6);grid on; 
axis([-Slen/2 Slen/2 -2 2]); 
title('Transmitted OFDM signal spectrum'); 
simulation_time = toc