您现在的位置:首页 >> 技术文章 >> MATLAB技术 >> 内容

MATLAB代做|FPGA代做--OFDM系统信道估计

时间:2018-10-29 14:44:54 点击:

  核心提示:OFDM系统信道估计...
echo off;                                     % 关闭回显 
clear all;                                    % 从内存中清除变量和函数
close all;                                    % 关闭所有图形
fprintf( '\n OFDM仿真\n \n') ;                % 设置显示格式
% --------------------------------------------- %
%                   参数定义                    %
% --------------------------------------------- %
IFFT_bin_length = 1024;                                % 发送端的IFFT变换长度, 接收端的FFT变换长度,R代表接受,T代表发送
carrier_count = 200;                                   % 子载波数
bits_per_symbol = 2;                                   % 位数/符号
symbols_per_carrier = 50;                              % 符号数/载波
cp_length = input('cp length = ');                     % 输入循环前缀长度   
d4 = input('d4 = ');                                   % 输入最大多径时延扩展
a4 = input('a4 = ');                                   % 输入最大多径时延扩展的系数
SNR = input('SNR = ');                                 % 输入信道信噪比(dB)                        
% --------------------------------------------- %
%                初始参数设置完毕               %
% --------------------------------------------- %
baseband_out_length = carrier_count * symbols_per_carrier * bits_per_symbol;        % 计算发送的二进制序列长度:基带传送长度=载波数× 符号数/载波×位数/符号
carriers = (1:carrier_count) + (floor(IFFT_bin_length/4) - floor(carrier_count/2));  % 载波(坐标):(1:200) + 256 - 100 = 157:356
conjugate_carriers = IFFT_bin_length - carriers + 2;                                 % 载波变换(坐标):1024 - (157:356) + 2 = 1026 - (157:356) = (869:670) 
%---------------------------------------
% 构造共轭时间-载波矩阵,以便应用所谓的RCC,Reduced Computational Complexity算法,即ifft之后结果为实数 
% 也可以用flipdim函数构造对称共轭矩阵
%---------------------------------------
% --------------------------------------------- %
%                    发送信号                   %
% --------------------------------------------- %
%---------------------------------------
% Generate a random binary output signal:
%   - a row of uniform random numbers (between 0 and 1), rounded to 0 or 1
%   - this will be the baseband signal which is to be transmitted.
%---------------------------------------
baseband_out = round(rand(1,baseband_out_length));
%---------------------------------------
% round:朝最近的整数取整,rand:产生均匀分布的随机数矩阵(1×baseband_out_length阶)               
% Convert to 'modulo N' integers where N = 2^bits_per_symbol
%   - this defines how many states each symbol can represent
%   - first, make a matrix with each column representing consecutive bits 
%     from the input stream and the number of bits in a column equal to the
%     number of bits per symbol
%   - then, for each column, multiply each row value by the power of 2 that 
%     it represents and add all the rows
%   - for example:  input 0 1 1 0 0 0 1 1 1 0
%                   bits_per_symbol = 2
%                   convert_matrix = 0 1 0 1 1
%                                    1 0 0 1 0
%                   modulo_baseband = 1 2 0 3 2
%---------------------------------------
convert_matrix = reshape(baseband_out,bits_per_symbol,length(baseband_out)/bits_per_symbol) ;
%---------------------------------------
% RESHAPE Change size. 把baseband_out变为M×N阶的矩阵
%    RESHAPE(X,M,N) returns the M-by-N matrix whose elements
%    are taken columnwise from X.  An error results if X does
%    not have M*N elements.
%---------------------------------------
for k = 1:(length(baseband_out)/bits_per_symbol)
  modulo_baseband(k) = 0 ;
for i = 1:bits_per_symbol
  modulo_baseband(k) =  modulo_baseband(k) + convert_matrix(i,k)*2^(bits_per_symbol - i) ;
end
end
%---------------------------------------
% Serial to Parallel Conversion       串并转换
%   - convert the serial modulo N stream into a matrix where each column 
%     represents a carrier and each row represents a symbol
%   - for example:
%      serial input stream = a b c d e f g h i j k l m n o p
%      parallel carrier distribution =
%                           C1/s1=a  C2/s1=b  C3/s1=c C4/s1=d
%                           C1/s2=e  C2/s2=f  C3/s2=g C4/s2=h
%                           C1/s3=i  C2/s3=j  C3/s3=k C4/s3=l
%                              .        .        .       .
%                              .        .        .       .
%---------------------------------------
carrier_matrix = reshape(modulo_baseband, carrier_count, symbols_per_carrier)';          % 生成时间-载波矩阵(carrier_count×symbols_per_carrier)
%---------------------------------------
% Apply differential coding to each carrier string
%   - append an arbitrary start symbol (let it be 0, that works for all 
%     values of bits_per_symbol) (note that this is done using a vertical
%     concatenation [x;y] of a row of zeros with the carrier matrix, sweet!)
%   - perform modulo N addition between symbol(n) and symbol(n-1) to get the 
%     coded value of symbol(n)
%   - for example:
%                   bits_per_symbol = 2 (modulo 4)
%                   symbol stream =  3 2 1 0 2 3
%                   start symbol = 0
%
%                   coded symbols = 0 + 3 = 3
%                                   3 + 2 = 11 = 1
%                                   1 + 1 = 2
%                                   2 + 0 = 2
%                                   2 + 2 = 10 = 0
%                                   0 + 3 = 3
%
%                   coded stream =  0 3 1 2 2 0 3
%---------------------------------------
% --------------------------------------------- %
%                   QDPSK调制                   %
% --------------------------------------------- %
carrier_matrix = [zeros(1,carrier_count);carrier_matrix];                                    % 添加一个差分调制的初始相位,为0
for i = 2:(symbols_per_carrier + 1)
    carrier_matrix(i,:) = rem(carrier_matrix(i,:)+carrier_matrix(i-1,:),2^bits_per_symbol);  % 差分调制(rem除后取余)
end
%---------------------------------------
% Convert the differential coding into a phase
%   - each phase represents a different state of the symbol
%   - for example:
%                 bits_per_symbol = 2 (modulo 4)
%                 symbols = 0 3 2 1
%                 phases =
%                           0 * 2pi/4 = 0 (0 degrees)
%                           3 * 2pi/4 = 3pi/2 (270 degrees)
%                           2 * 2pi/4 = pi (180 degrees)
%                           1 * 2pi/4 = pi/2 (90 degrees)
%---------------------------------------
carrier_matrix = carrier_matrix * ((2*pi)/(2^bits_per_symbol));                              % 产生差分相位
%---------------------------------------
% Convert the phase to a complex number
%   - each symbol is given a magnitude of 1 to go along with its phase 
%     (via the ones(r,c) function)
%   - it is then converted from polar to cartesian (complex) form
%   - the result is 2 matrices, X with the real values and Y with the imaginary
%   - each X column has all the real values for a carrier, and each Y column 
%     has the imaginary values for a carrier
%   - a single complex matrix is then generated taking X for real and 
%     Y for imaginary
%---------------------------------------
[X,Y] = pol2cart(carrier_matrix, ones(size(carrier_matrix,1),size(carrier_matrix,2)));       % 由极坐标向复数坐标转化 第一参数为相位 第二参数为幅度
%---------------------------------------
% Carrier_matrix contains all the phase information and all the amplitudes are the same,‘1’     
% 函数说明:极或柱坐标变为直角坐标
%  POL2CART Transform polar to Cartesian coordinates.
%    [X,Y] = POL2CART(TH,R) transforms corresponding elements of data
%    stored in polar coordinates (angle TH, radius R) to Cartesian
%    coordinates X,Y.  The arrays TH and R must the same size (or
%    either can be scalar).  TH must be in radians.
% [X,Y,Z] = POL2CART(TH,R,Z) transforms corresponding elements of
%    data stored in cylindrical coordinates (angle TH, radius R, height Z) 
%    to Cartesian coordinates X,Y,Z. The arrays TH, R, and Z must be
%    the same size (or any of them can be scalar).  TH must be in radians.
%---------------------------------------
complex_carrier_matrix = complex(X,Y);                     
%---------------------------------------
% 函数说明:
%  COMPLEX Construct complex result from real and imaginary parts.
%    C = COMPLEX(A,B) returns the complex result A + Bi, where A and B are
%    identically sized real N-D arrays, matrices, or scalars of the same data type.
%
% Assign each carrier to its IFFT bin
%   - each row of complex_carrier_matrix represents one symbol period, with
%     a symbol for each carrier
%   - a matrix is generated to represent all IFFT bins (columns) and all 
%     symbols (rows)
%   - the phase modulation for each carrier is then assigned to the 
%     appropriate bin
%   - the conjugate of the phase modulation is then assigned to the 
%     appropriate bin
%      - the phase modulation bins and their conjugates are symmetric about 
%        the Nyquist frequency in the IFFT bins
%      - since the first bin is DC, the Nyquist Frequency is located 
%        at (number of bins/2) + 1
%      - symmetric conjugates are generated so that when the signal is 
%        transformed to the time domain, the time signal will be real-valued
%   - example
%      - 1024 IFFT bins
%      - bin 513 is the center (symmetry point)
%      - bin 1 is DC
%      - bin 514 is the complex conjugate of bin 512
%      - bin 515 is the complex conjugate of bin 511
%      - ....
%      - bin 1024 is the complex conjugate of bin 2 (if all bins 
%        were used as carriers)
%      - So, bins 2-512 map to bins 1024-514
%---------------------------------------
%-------------------------------------- %
%                               添加训练序列                                  %
%-------------------------------------- % 
training_symbols = [ 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 ...
-j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 ...
1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 ...
-1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j ...
-1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 1 j j 1 -1 -j -j -1 ]; 
%---------------------------------------
% 25 times "1 j j 1" 
% 25 times "-1 -j -j -1"
% totally 200 symbols as a row
%---------------------------------------
training_symbols = cat(1, training_symbols, training_symbols) ;  
training_symbols = cat(1, training_symbols, training_symbols) ;   
%---------------------------------------
% Production of 4 rows of training_symbols
% 函数说明:串接成高维数组
% CAT Concatenate arrays.
%     CAT(DIM,A,B) concatenates the arrays A and B along the dimension DIM.  
%     CAT(2,A,B) is the same as [A,B].
%     CAT(1,A,B) is the same as [A;B].
%     B = CAT(DIM,A1,A2,A3,A4,...) concatenates the input
%     arrays A1, A2, etc. along the dimension DIM.
%---------------------------------------
complex_carrier_matrix = cat(1, training_symbols, complex_carrier_matrix);
%---------------------------------------
% 训练序列与数据合并,串接成高维数组
% block-type pilot symbols
%---------------------------------------
IFFT_modulation = zeros(4 + symbols_per_carrier + 1, IFFT_bin_length);
%---------------------------------------
% Here a row vector of zeros is between training symbols and data symbols!!! 
% 4 training symbols and 1 zero symbol
% every OFDM symbol takes a row of "IFFT_modulation" 
%---------------------------------------
IFFT_modulation(:,carriers) = complex_carrier_matrix;
IFFT_modulation(:,conjugate_carriers) = conj(complex_carrier_matrix);
%---------------------------------------             
% Find the indices of zeros
% Description
% ZC = conj(Z) returns the complex conjugate of the elements of Z. 
% AlgorithmIf Z is a complex array: conj(Z) = real(Z) - i*imag(Z)
%---------------------------------------
time_wave_matrix = ifft(IFFT_modulation');                                     % 进行IFFT操作
time_wave_matrix = time_wave_matrix';                                          % 进行矩阵转置
cp_add = zeros(4 + symbols_per_carrier + 1,cp_length);
 for i = 1:4 + symbols_per_carrier + 1
    cp_add(i,:) = time_wave_matrix(i,(IFFT_bin_length - cp_length + 1) : IFFT_bin_length);
end
time_wave_matrix_cp = [cp_add,time_wave_matrix];
for i = 1: 4 + symbols_per_carrier + 1                     % windowed_time_wave_matrix(i,:) = real(time_wave_matrix(i,:)) .* hamming(IFFT_bin_length)';
    windowed_time_wave_matrix_cp(i,:) = real(time_wave_matrix_cp(i,:));
end
%---------------------------------------
% Serialize the modulating waveform
%   - sequentially take each row of windowed_time_wave_matrix and construct a row vector
%   - the row vector will be the modulating signal
%   - note that windowed_time_wave_matrix is transposed, this is to account for the way the
%     Matlab 'reshape' function works (reshape takes the columns of the target matrix and 
%     appends them sequentially) 
% get the real part of the result of IFFT
% 这一步可以省略,因为IFFT结果都是实数
% 由此可以看出,只是取了IFFT之后载波上的点,并未进行CP的复制和添加end
%---------------------------------------
ofdm_modulation = reshape(windowed_time_wave_matrix_cp', 1, (IFFT_bin_length + cp_length)*(4 + symbols_per_carrier+1));   % P2S operation
%---------------------------------------
Tx_data = ofdm_modulation;
%---------------------------------------
%                                  信道模拟                                    
%---------------------------------------
% The channel model is Gaussian (AWGN) only 
%   - Rayleigh fading would be a useful addition
%---------------------------------------
d1 = 40; a1 = 0.2; d2 = 50; a2 = 0.3; d3 = 60; a3 = 0.4; 
%d4 = 160; a4 = 0.9;
copy1 = zeros(size(Tx_data)) ;
for i = 1 + d1: length(Tx_data)
copy1(i) = a1*Tx_data( i - d1) ;
end
copy2 = zeros(size(Tx_data) ) ;
for i = 1 + d2: length( Tx_data)
copy2(i) = a2*Tx_data( i - d2) ;
end
copy3 = zeros(size(Tx_data) ) ;
for i = 1 + d3: length(Tx_data)
copy3(i) = a3*Tx_data ( i - d3) ;
end
copy4 = zeros(size(Tx_data) ) ;
for i = 1 + d4: length( Tx_data)
copy4(i) = a4*Tx_data(i - d4) ;
end
Tx_data = Tx_data + copy1 + copy2 + copy3 + copy4;           % 4 multi-paths
Tx_signal_power = var(Tx_data);
%-------------------------------------------------------------------------
% 函数说明:
%  VAR Variance.
%     For vectors, Y = VAR(X) returns the variance of the values in X.  For
%     matrices, Y is a row vector containing the variance of each column of
%     X.
%-------------------------------------------------------------------------
linear_SNR = 10^(SNR/10);
noise_sigma = Tx_signal_power/linear_SNR;
noise_scale_factor = sqrt(noise_sigma);
%
noise = randn(1, length(Tx_data))*noise_scale_factor;
Rx_Data = Tx_data + noise;
%Rx_Data = Tx_data;
%-------------------------------------------------------------------------
% 函数说明:产生正态分布的随机函数
% Y = randn(m,n) or Y = randn([m n]) returns an m-by-n matrix of random
% entries.
% The randn function generates arrays of random numbers whose elements are
% normally distributed with mean 0 and variance 1.
%-------------------------------------------------------------------------
% --------------------------------------------- %
%                  信号接收                      %
% --------------------------------------------- %
%-------------------------------------------------------------------------
Rx_Data_matrix_cp = reshape(Rx_Data, IFFT_bin_length + cp_length, 4 + symbols_per_carrier + 1);
Rx_Data_matrix_cp = Rx_Data_matrix_cp';
Rx_Data_matrix = zeros(4 + symbols_per_carrier + 1,IFFT_bin_length);
for i = 1:4 + symbols_per_carrier + 1
    Rx_Data_matrix(i,:) = Rx_Data_matrix_cp(i,(cp_length + 1):(IFFT_bin_length + cp_length));
end
Rx_Data_matrix = Rx_Data_matrix';
%-------------------------------------------------------------------------
% Transform each symbol from time to frequency domain
%   - take the fft of each column
%-------------------------------------------------------------------------
Rx_spectrum = fft(Rx_Data_matrix);
%-------------------------------------------------------------------------
%  Suppose precise synchronazition between Tx and Rx
Rx_carriers = Rx_spectrum(carriers,:)';
Rx_training_symbols = Rx_carriers( (1: 4) , : ) ;
Rx_carriers = Rx_carriers((5: 55), : ) ;
% --------------------------------------------- %
%                    信道估计                    %
% --------------------------------------------- %
Rx_training_symbols = Rx_training_symbols./ training_symbols;
Rx_training_symbols_deno = Rx_training_symbols.^2;
Rx_training_symbols_deno = Rx_training_symbols_deno(1,:)+Rx_training_symbols_deno(2,:)+Rx_training_symbols_deno(3,:)+Rx_training_symbols_deno(4,:) ;
Rx_training_symbols_nume = Rx_training_symbols(1, : ) +Rx_training_symbols(2, : ) + Rx_training_symbols(3, : ) +Rx_training_symbols(4, : ) ;
Rx_training_symbols_nume = conj(Rx_training_symbols_nume) ;
%-------------------------------------------------------------------------
% 取4个向量的导频符号是为了进行平均优化
% 都是针对 “行向量”即单个的OFDM符号 进行操作
% 原理:寻求1/H,对FFT之后的数据进行频域补偿
% 1/H = conj(H)/H^2 because H^2 = H * conj(H) 
%-------------------------------------------------------------------------
Rx_training_symbols = Rx_training_symbols_nume./Rx_training_symbols_deno;
Rx_training_symbols_2 = cat(1, Rx_training_symbols,Rx_training_symbols) ;
Rx_training_symbols_4 = cat(1, Rx_training_symbols_2,Rx_training_symbols_2) ;
Rx_training_symbols_8 = cat(1, Rx_training_symbols_4,Rx_training_symbols_4) ;
Rx_training_symbols_16 = cat(1, Rx_training_symbols_8, Rx_training_symbols_8) ;
Rx_training_symbols_32 = cat(1, Rx_training_symbols_16, Rx_training_symbols_16) ;
Rx_training_symbols_48 = cat(1, Rx_training_symbols_32, Rx_training_symbols_16) ;
Rx_training_symbols_50 = cat(1, Rx_training_symbols_48, Rx_training_symbols_2) ;
Rx_training_symbols = cat(1, Rx_training_symbols_50,Rx_training_symbols) ;
Rx_carriers = Rx_training_symbols.*Rx_carriers;
%-------------------------------------------------------------------------
% 进行频域单抽头均衡
%-------------------------------------------------------------------------
Rx_phase = angle(Rx_carriers)*(180/pi);
phase_negative = find(Rx_phase < 0);
%-------------------------------------------------------------------------
% 函数说明:找出非零元素的索引号
%     FIND   Find indices of nonzero elements.
%     I = FIND(X) returns the linear indices of the array X that are nonzero.
%     X may be a logical expression. Use IND2SUB(I,SIZE(X)) to calculate
%     multiple subscripts from the linear indices I.
%---------------------------------------------------------------------------
Rx_phase( phase_negative );
Rx_phase(phase_negative) = rem(Rx_phase(phase_negative) + 360, 360) ;
% 把负的相位转化为正的相位
Rx_decoded_phase = diff(Rx_phase) ;
%  这也是为什么要在前面加上初始相位的原因 
% “Here a row vector of zeros is between training symbols and data symbols!!!”
%-------------------------------------------------------------------------
% 函数说明:
%     DIFF Difference and approximate derivative.
%     DIFF(X), for a vector X, is [X(2)-X(1)  X(3)-X(2) ... X(n)-X(n-1)].
%     DIFF(X), for a matrix X, is the matrix of row differences,
%        [X(2:n,:) - X(1:n-1,:)].
%------------------------Test Codes --------------------------------------
%  a = [1 2 3; 4 5 6; 7 8 9; 10 11 12];
%  b = a;
%      for i = 2:4
%        b(i,:) = b(i-1,:) + b(i,:);
%      end
%  c = diff(b);
%-----------------------Test Result --------------------------------------
%  a =                  Modulating signal
%      1     2     3
%      4     5     6
%      7     8     9
%     10    11    12
%  b =                  Modulated signal
%      1     2     3
%      5     7     9
%     12    15    18
%     22    26    30
%  c =                  Recovered signal
%      4     5     6
%      7     8     9
%     10    11    12
% ----------------------------------------------------------------------------
%   Name                   Size                    Bytes  Class
%   Rx_phase              51x200                   81600  double array
%   Rx_decoded_phase      50x200                   80000  double array
% ---------------------------------------------------------------------------- 
phase_negative = find(Rx_decoded_phase < 0) ;
Rx_decoded_phase(phase_negative)= rem(Rx_decoded_phase(phase_negative) + 360, 360) ;
% ---------------------------------------------------------------------------- 
% Extract phase differences (from the differential encoding)
%   - the matlab diff( ) function is perfect for this operation
%   - again, normalize the result to be between 0 and 359 degrees
% 再次把负的相位转化为正的相位
% ---------------------------------------------------------------------------- 
% --------------------------------------------- %
%                   QDPSK解调                   %
% --------------------------------------------- %
base_phase = 360 /2^bits_per_symbol;
delta_phase = base_phase /2;
Rx_decoded_symbols = zeros(size(Rx_decoded_phase,1),size(Rx_decoded_phase,2)) ;
%
for i = 1: (2^bits_per_symbol - 1)
center_phase = base_phase*i;
plus_delta = center_phase + delta_phase;  % Decision threshold 1
minus_delta = center_phase - delta_phase; % Decision threshold 2
decoded = find((Rx_decoded_phase <= plus_delta)&(Rx_decoded_phase > minus_delta)) ;
Rx_decoded_symbols(decoded) = i;
end
% ---------------------------------------------------------------------------- 
%  仅仅对三个区域进行判决
%  剩下的区域就是零相位的空间了
%  这个区域在定义解调矩阵时已经定义为零
% ---------------------------------------------------------------------------- 
Rx_serial_symbols = reshape(Rx_decoded_symbols',1,size(Rx_decoded_symbols, 1)*size(Rx_decoded_symbols,2)) ;
for i = bits_per_symbol: -1: 1
if i ~= 1
Rx_binary_matrix(i, : ) = rem(Rx_serial_symbols, 2) ;
Rx_serial_symbols = floor(Rx_serial_symbols/2) ;
else
Rx_binary_matrix( i, : ) = Rx_serial_symbols;
end
end                                                             % Integer to binary
baseband_in = reshape(Rx_binary_matrix, 1,size(Rx_binary_matrix, 1)*size(Rx_binary_matrix, 2) ) ;
% --------------------------------------------- %
%                   误码率计算                   %
% --------------------------------------------- %
bit_errors = find(baseband_in ~= baseband_out) ;            % find的结果 其每个元素为满足逻辑条件的输入向量的标号,其向量长度也就是收发不一样的bit的个数
bit_error_count = size(bit_errors, 2) ;
total_bits = size( baseband_out, 2) ;
bit_error_rate = bit_error_count/ total_bits;

fprintf ( '%f \n',bit_error_rate) ;


联系:highspeedlogic

QQ :1224848052

微信:HuangL1121

邮箱:1224848052@qq.com

网站:http://www.mat7lab.com/

网站:http://www.hslogic.com/

微信扫一扫:

作者:OFDM系统信道估计 来源:OFDM系统信道估计
  • 您是如何找到本站的?
  • 百度搜索
  • Google搜索
  • 查阅资料过程中
  • 论坛发现
  • 百度贴吧发现
  • 朋友介绍
本站最新成功开发工程项目案例
相关文章
  • 没有相关文章
相关评论
发表我的评论
  • 大名:
  • 内容:
  • matlab代做|matlab专业代做|matlab淘宝代做(www.hslogic.com) © 2018 版权所有 All Rights Reserved.
  • Email:highspeed_logic@163.com 站长QQ: 1224848052

    专业代做/代写/承接、MATLAB、SIMULINK、FPGA项目、博士/硕士/本科毕业设计、课题设计、论文,毕业论文,Coursework、Eassy、Assignment