OFDM信道估计MATLAB代码

1、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(&

2、#39;cp length = ' % 输入循环前缀长度 d4 = input('d4 = ' % 输入最大多径时延扩展 a4 = input('a4 = ' % 输入最大多径时延扩展的系数 SNR = input('SNR = ' % 输入信道信噪比(dB % - % % 初始参数设置完毕 % % - % baseband_out_length = carrier_count * symbols_per_carrier * bits_per_symbol; % 计算发送的二进制序列长度：基带传送长度=载波数× 符号数/载波&#

3、215;位数/符号 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算

4、法，即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; % % ro

5、und：朝最近的整数取整，rand:产生均匀分布的随机数矩阵（1×baseband_out_length阶） % Convert to 'modulo N' integers where N = 2bits_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

6、 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

7、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

8、 = 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 % repre

9、sents 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 = res

10、hape(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 (not e that this is done using a ver

11、tical % 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 =

12、 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(

13、i-1,:,2bits_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

14、/4 = pi (180 degrees % 1 * 2pi/4 = pi/2 (90 degrees % carrier_matrix = carrier_matrix * (2*pi/(2bits_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

15、 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 imagi

16、nary % 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(

17、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 c

18、ylindrical 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 pa

19、rts. % 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 %

20、 - 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 conjuga

21、tes 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 % - exam

22、ple % - 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 10

23、24-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

24、-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 time

25、s "-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 th

26、e 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; % % 训练序列与数据合并，串接成高维数

27、组 % 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(:,

28、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 = iff

29、t(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_wa

30、ve_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_wav

31、e_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 % g

32、et 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 mode

33、l 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(

34、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 ; f or i = 1 + d4: length( Tx_data copy4(i = a4*Tx_data(i - d4 ; end Tx_data = Tx_data + copy1 + copy2 + copy3 + copy4; % 4 multi-path

35、s 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 = sqr

36、t(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 % normal

37、ly 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_pe

38、r_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 synchronazit

39、ion 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

40、_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_traini

41、ng_symbols_nume ; %- - % 取4个向量的导频符号是为了进行平均优化 % 都是针对 “行向量”即单个的OFDM符号 进行操作 % 原理：寻求1/H，对FFT之后的数据进行频域补偿 % 1/H = conj(H/H2 because H2 = 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_train

42、ing_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

43、_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; %- % 进行频域单抽

44、头均衡 %- 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 sub

45、scripts 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 D

46、ifference 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;

47、 %-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 = fi