1、echo off; % 關(guān)閉回顯 clear all; % 從內(nèi)存中清除變量和函數(shù)close all; % 關(guān)閉所有圖形fprintf( n OFDM仿真n n) ; % 設(shè)置顯示格式% - % 參數(shù)定義 % - %IFFT_bin_length = 1024; % 發(fā)送端的IFFT變換長度, 接收端的FFT變換長度，R代表接受，T代表發(fā)送carrier_count = 200; % 子載波數(shù)bits_per_symbol = 2; % 位數(shù)/符號symbols_per_carrier = 50; % 符號數(shù)/載波cp_length = input(cp length = ); % 輸入循環(huán)前

2、綴長度 d4 = input(d4 = ); % 輸入最大多徑時延擴展a4 = input(a4 = ); % 輸入最大多徑時延擴展的系數(shù)SNR = input(SNR = ); % 輸入信道信噪比(dB) % - % 初始參數(shù)設(shè)置完畢 % - %baseband_out_length = carrier_count * symbols_per_carrier * bits_per_symbol; % 計算發(fā)送的二進(jìn)制序列長度：基帶傳送長度=載波數(shù) 符號數(shù)/載波位數(shù)/符號carriers = (1:carrier_count) + (floor(IFFT_bin_length/4) - flo

3、or(carrier_count/2); % 載波（坐標(biāo)）：(1:200) + 256 - 100 = 157:356conjugate_carriers = IFFT_bin_length - carriers + 2; % 載波變換（坐標(biāo)）：1024 - (157:356) + 2 = 1026 - (157:356) = (869:670) % 構(gòu)造共軛時間載波矩陣，以便應(yīng)用所謂的RCC，Reduced Computational Complexity算法，即ifft之后結(jié)果為實數(shù) % 也可以用flipdim函數(shù)構(gòu)造對稱共軛矩陣% - % 發(fā)送信號 % - % Generate a ra

4、ndom 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：朝最近的整數(shù)取整，rand:產(chǎn)生均勻分布的隨機數(shù)矩陣（1baseband_out_length階） % Convert to modulo N integers

5、 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 bits in a column equal to the% number of bits per symbol% - then, for each column, multiply each row

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

7、Change size. 把baseband_out變?yōu)镸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_

8、baseband(k) + convert_matrix(i,k)*2(bits_per_symbol - i) ;endend% Serial to Parallel Conversion 串并轉(zhuǎn)換% - 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% p

9、arallel 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_countsymbols_per_carrier）% Apply differential coding t

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

11、 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調(diào)制 % - %carrier_matrix = zeros(1,carrier_

12、count);carrier_matrix; % 添加一個差分調(diào)制的初始相位，為0for i = 2:(symbols_per_carrier + 1) carrier_matrix(i,:) = rem(carrier_matrix(i,:)+carrier_matrix(i-1,:),2bits_per_symbol); % 差分調(diào)制（rem除后取余）end% Convert the differential coding into a phase% - each phase represents a different state of the symbol% - for example

13、:% 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)/(2bits_per_symbol); % 產(chǎn)生差分相位% Convert the phase to a complex number% - each

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

15、olumn % 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); % 由極坐標(biāo)向復(fù)數(shù)坐標(biāo)轉(zhuǎn)化 第一參數(shù)為相位 第二參數(shù)為幅度% Carrier_matrix contains all the phase information and

16、all the amplitudes are the same,1 % 函數(shù)說明：極或柱坐標(biāo)變?yōu)橹苯亲鴺?biāo)% 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

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

18、 radians.%complex_carrier_matrix = complex(X,Y); % 函數(shù)說明：% 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

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

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

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

22、the complex conjugate of bin 2 (if all bins % were used as carriers)% - So, bins 2-512 map to bins 1024-514% % 添加訓(xùn)練序列 % % 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 -

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

24、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% 函數(shù)說明：串接成高

25、維數(shù)組% 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, traini

26、ng_symbols, complex_carrier_matrix);% 訓(xùn)練序列與數(shù)據(jù)合并，串接成高維數(shù)組% 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

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

28、) - i*imag(Z)%time_wave_matrix = ifft(IFFT_modulation); % 進(jìn)行IFFT操作time_wave_matrix = time_wave_matrix; % 進(jìn)行矩陣轉(zhuǎn)置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);endtime_wav

29、e_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 r

30、ow 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 sequen

31、tially) % get the real part of the result of IFFT% 這一步可以省略，因為IFFT結(jié)果都是實數(shù)% 由此可以看出，只是取了IFFT之后載波上的點，并未進(jìn)行CP的復(fù)制和添加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 i

32、s 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) ;endcopy2 = zeros(size(Tx_data) ) ;for i = 1 + d2: length( Tx_

33、data)copy2(i) = a2*Tx_data( i - d2) ;endcopy3 = zeros(size(Tx_data) ) ;for i = 1 + d3: length(Tx_data)copy3(i) = a3*Tx_data ( i - d3) ;endcopy4 = zeros(size(Tx_data) ) ;for i = 1 + d4: length( Tx_data)copy4(i) = a4*Tx_data(i - d4) ;endTx_data = Tx_data + copy1 + copy2 + copy3 + copy4; % 4 multi-path

34、sTx_signal_power = var(Tx_data);%-% 函數(shù)說明：% 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_

35、sigma);%noise = randn(1, length(Tx_data)*noise_scale_factor;Rx_Data = Tx_data + noise;%Rx_Data = Tx_data;%-% 函數(shù)說明：產(chǎn)生正態(tài)分布的隨機函數(shù)% 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

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

37、(i,:) = Rx_Data_matrix_cp(i,(cp_length + 1):(IFFT_bin_length + cp_length);endRx_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 RxRx_carriers = Rx

38、_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_trai

39、ning_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個向量的導(dǎo)頻符號是為了

40、進(jìn)行平均優(yōu)化% 都是針對 “行向量”即單個的OFDM符號 進(jìn)行操作% 原理：尋求1/H，對FFT之后的數(shù)據(jù)進(jìn)行頻域補償% 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_training_symbols_4 = cat(1, Rx_training_symbol

41、s_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_5