完整的OFDM系统的仿真实现
%窗函数子程序，子程序名称：recoswindow.m
function[rcosw]=rcoswindow(beta,Ts)
%输入参数：beta为升余弦窗关键系数，Ts为IFFT长度家循环前缀长度t=0:(1+beta)*Ts;
rcosw=zeros(1,(1+beta)*Ts);
%计算升余弦窗，共有三部分
for i=1:beta*Ts;
rcosw(i)=0.5+0.5*cos(pi+t(i)*pi/(beta*Ts));%计算升余弦窗的第一部分end
rcosw(beta*Ts+1:Ts)=1;%计算升余弦窗低二部分
for j=Ts:(1+beta)*Ts+1;
rcosw(j-1)=0.5+0.5*cos((t(j)-Ts)*pi/(beta*Ts));
%计算升余弦窗第三部分
end
rcosw=rcosw';%转换为列矢量
%将16QAM信号的解调子程序，子程序的名称：demoduqam16.m
end
%16QAM调制子程序，子程序名称为qam16.m
%将二进制数目流转换为16QAM信号
function [complex_qam_data]=qam16(bitdata)
%输入参数：bitdata为二进制数码流
%输出参数：complex_qam_data为16QAM副信号
X1=reshape(bitdata,4,length(bitdata)/4)';%将二进制数码流以4bitte分段
d=1;
%转换4bit二进制码为十进制码1~16，生态农场mapping映射表中的索引for i=1:length(bitdata)/4;
for j=1:4
X1(i,j)=X1(i,j)*(2^(4-j));
end
source(i,1)=1+sum(X1(i,:));
end
%16QAM映射表，改表中存放的16对，没对两个实数，表示星座位置
mapping=[-3*d 3*d;-d 3*d;d 3*d;3*d 3*d;-3*d d;-d d;d d;3*d d;-3*d -d;-d -d;d -d;3*d -d;-3*d -3*d;-d -3*d;d -3*d;3*d -3*d];
for i=1:length(bitdata)/4
qam_data(i,:)=mapping(source(i),:);%数据映射
end
complex_qam_data=complex(qam_data(:,1),qam_data(:,2));
%组合为负数形式，形成16QAM信号
end
%将16QAM信号的解调子程序，子程序的名称：demoduqam16.m
%该子程序测试
function [demodu_bit_symble]=demoduqam16(Rx_serial_complex_symbols)
%输入参数为：Rx_serial_complex_symbols为接收端接收到的复16QAM信号%输出参数：demodu_bit_symble为二进制数码流
complex_symbols=reshape(Rx_serial_complex_symbols,length(Rx_serial_comple x_symbols),1);
d=1;
mapping=[-3*d 3*d;-d 3*d;d 3*d;3*d 3*d;-3*d d;-d d;d d;3*d d;-3*d -d;-d -d;d -d;3*d -d;-3*d -3*d;-d -3*d;d -3*d;3*d -3*d];
complex_mapping=complex(mapping(:,1),mapping(:,2));
%将数据映射表中转换为16QAM信号，即3组合为复数；
for i=1:length(Rx_serial_complex_symbols);
for j=1:16;
metrics(j)=abs(complex_symbols(i,1)-complex_mapping(j,1));
end
[min_metricdecode_symble(i)]=min(metrics);
%将接收数据与标准16QAM信号比，找到差最小的，将其对应恢复成标准的16QAM信号
end
decode_bit_symble=de2bi((decode_symble-1)','left-msb');
%将16QAM转为二进制
demodu_bit_symble=reshape(decode_bit_symble',1,length(Rx_serial_complex_sy mbols)*4);%转换为一行
end
baseband_out_length=16000;
rand('twister',0);
baseband_out=round(rand(1,baseband_out_length));
%产生16000bit待传输的二进制比特流。这里存放的是发送的二进制信号与后面解调后的二进制信号比较，可以计算误码率。
%16QAM调制病绘制星座图。
complex_carrier_matrix=qam16(baseband_out);
figure(1);
plot(complex_carrier_matrix,'*r');%绘制16QAM星座图
title('16QAM调制后星座图');
grid on ;
%16QAM调制子程序，子程序名称为qam16.m
%将二进制数目流转换为16QAM信号
%QAM16测试
%输入参数：bitdata为二进制数码流
%输出参数：complex_qam_data为16QAM复信号
close all;
clear all;
baseband_out_length=16000;
rand('twister',0);
bitdata=round(rand(1,baseband_out_length));
complex_carrier_matrix=qam16(bitdata);
% X1=reshape(bitdata,4,length(bitdata)/4)';%将二进制数码流以4bite分段
% d=1;
% %转换4bit二进制码为十进制码1~16，生成mapping映射表中的索引
% for i=1:length(bitdata)/4;
% for j=1:4
% X1(i,j)=X1(i,j)*(2^(4-j));
% end
% source(i,1)=1+sum(X1(i,:));
% end
% %16QAM映射表,该表中存放的16对，没对两个实数，表示星座位置
% mapping=[-3*d 3*d;-d 3*d;d 3*d;3*d 3*d;-3*d d;-d d;d d;3*d d;-3*d -d;-d -d;d -d;3*d -d;-3*d -3*d;-d -3*d;d -3*d;3*d -3*d];
% for i=1:length(bitdata)/4
% qam_data(i,:)=mapping(source(i),:);%数据映射
% end
% complex_qam_data=complex(qam_data(:,1),qam_data(:,2)); % %组合为负数形式，形成16QAM信号
% figure(1);
% plot(complex_qam_data,'*r');%绘制16QAM星座图
% title('16QAM调制后星座图');
% grid on ;
figure(1);
plot(complex_carrier_matrix,'*r');%绘制16QAM星座图
title('16QAM调制后星座图');
grid on ;
% 程序功能
% （1）通过对OFDM信号各个子载波赋共轭对称的数据产生一个实OFDM 符号。
% （2）给OFDM符号加循环前缀与循环后缀
% （3）给OFDM符号加窗。在程序中加入的是升余弦窗，可以通过改变升余弦窗的滚降系数，观察加入不同升余弦窗对OFDM信号频谱的的影响% （4）信道采用嘉兴高斯白噪声信道。可以通过改变信噪比来改变信道环境，从而在接受端通过误码率或是星座图观察信道对OFDM信号传输的影响% （5）去除循环前缀与循环后缀，对OFDM信号进行解调
% 程序代码：
clear all;
close all;
carrier_count=200;
%这个程序中OFDM子载波个数为512，其中400即carrier_count *2为数据符号，其余赋0值
symbols_per_carrier=20;
%每个子载波上的符号数，在这里即是OFDM符号的个数
bits_per_symbol=4;
%OFDM符号的每个子载波上传输的比特数，4bit通常采用的16QAM调制IFFT_bin_length=512;
%FFT的长度，也即一个OFDM符号的子载波个数；
PrefixRatio=1/4;
%循环前缀的比值，循环前缀与OFDM符号长度的壁纸，通常在1/6~1/4之间GI=PrefixRatio * (IFFT_bin_length);%保护间隔的长度，这里为128
beta=1/32;%升余弦窗的滚降系数
GIP=beta*(IFFT_bin_length+GI);%循环后缀的长度，这里为20；
SNR=16;%本程序考虑加性高斯白噪声信道，这里信噪比为30db
%===============================OFDM信号产生=====================
baseband_out_length=carrier_count*symbols_per_carrier * bits_per_symbol;
%计算传输数据总的比特数，为200*20*4
%bits=16000bits。16000bits构成20个OFDM符号，每个OFDM符号200个子载波，每个子载波传输4bit信息
carriers=(1:carrier_count)+(floor(IFFT_bin_length/4)-floor(carrier_count/2));
%计算OFDM符号的在载波的序号，carriers 中存放的序号是29~228；
conjugate_carriers=IFFT_bin_length-carriers +2;
% conjugate_carriers=IFFT_bin_length-carriers -2;
%====================================================这个序号有问题
%计算OFDM符号子载波的序号，conjugate_carriers中存放的序号是282~481；
rand('twister',0);
baseband_out=round(rand(1,baseband_out_length));
%产生16000bit待传输的二进制比特流。这里存放的是发送的二进制信号与后面解调后的二进制信号比较，可以计算误码率。
%16QAM调制病绘制星座图。
complex_carrier_matrix=qam16(baseband_out);
%调用子程序qam16进行16QAM调制。将base_band中的二进制比特流，每4bit转换为一个16QAM信号，即将二进制比特流每4bit转换为-3-3j，-3+3j,3-3j,3+3j, %-1-3j,-1+3j,1-3j,1+3j,-3-j,-3+j,3-j,3+j,-1-j,-1+j,1-j,1+j中的一个。转换后complex_carrier_matrix为1*4000矩阵。
complex_carrier_matrix=reshape(complex_carrier_matrix',carrier_count,symbols_ per_carrier);
%转换complex_carrier_matrix中的数据为carrier_coun6t*symbols_per_carrier 矩阵，这里为20*200矩阵。
figure(1);
plot(complex_carrier_matrix,'*r');%绘制16QAM星座图
title('16QAM调制后星座图');
grid on ;
%IFFT,即进行OFDM调制
IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);
%将symbols_per_carrier*IFFT_bin_length矩阵赋0值，这里讲20*512矩阵赋0值。这里512是IFFT的长度，也是OFDM符号的子载波个数
%这里明显有问题,以下矫正
% IFFT_modulation(:,carriers)=complex_carrier_matrix;
IFFT_modulation(:,carriers)=complex_carrier_matrix';
%将20*200的complex_carrier_matrix的数据赋给IFFT_modulation的第29~228列，即给512个子载波中的第29~229个子载波赋值
% IFFT_modulation(:,conjugate_carriers)=conj(complex_carrier_matrix);
IFFT_modulation(:,conjugate_carriers)=conj(complex_carrier_matrix');
%将20*200的complex_carrier_matrix的数据付给512个子载波中的第282~481个子载波，这段程序构造了512个在载波的OFDM符号，并且各个%子载波上的数据是共轭对称的。这样做的目的是经过IFFT后形成的OFDM 符号均为实数。另外，在512个子载波中，仅有400个子载波为数据，其余为%0值，相当于补零。补零的目的是通常IFFT的长度应该为2的证书次幂。
signal_after_IFFT=ifft(IFFT_modulation,IFFT_bin_length,2);
%IFFT实现OFDM调制
time_wave_matrix=signal_after_IFFT;
figure(2);
plot(0:IFFT_bin_length-1,time_wave_matrix(2,:));
%画一个OFDM信号的时域表现
axis([0,512,-0.4,0.4]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('OFDM Time Signal,One Symbol Period');
%添加循环前缀与循环后缀
XX=zeros(symbols_per_carrier,IFFT_bin_length +GI+GIP);
%IFFT_bin_length+GI+GIP为OFDM，循环前缀，循环后缀长度之和；for k=1:symbols_per_carrier;
for i=1:IFFT_bin_length;
XX(k,i+GI)=signal_after_IFFT(k,i);
end
for i=1:GI
XX(k,i)=signal_after_IFFT(k,i+IFFT_bin_length-GI);
%添加循环前缀
end
for j=1:GIP
XX(k,IFFT_bin_length+GI+j)=signal_after_IFFT(k,j);
%添加循环后缀
end
end
time_wave_matrix_cp=XX;%带循环前缀和循环后缀的OFDM符号figure(3)
plot(0:length(time_wave_matrix_cp)-1,time_wave_matrix_cp(2,:));
%画带循环前缀与循环后缀的OFDM信号的时域波形、
axis([0,600,-0.3,0.3]);
grid on ;
ylabel('Ampelitude');
xlabel('Time');
title('OFDM Time Signal with CP,One Symbol Period');
%OFDM符号加窗
windowed_time_wave_matrix_cp=zeros(1,IFFT_bin_length+GI+GIP);
for i=1:symbols_per_carrier
windowed_time_wave_matrix_cp(i,:)=real(time_wave_matrix_cp(i,:)).*rcoswindow(bet a,IFFT_bin_length+GI)';
%调用rcoswindow产生升余弦窗，对待循环前缀与循环后缀的OFDM符号加窗
end
figure(4);
plot(0:IFFT_bin_length-1+GI+GIP,windowed_time_wave_matrix_cp(2,:));
%画加窗后的OFDM符号
axis([0,700,-0.2,0.2]);
grid on ;
ylabel('Amplitude');
xlabel('Time');
title('OFDM Time Signal Apply a Window,One Symbol Period');
%生成发送信号，并串转换
windowed_Tx_data=zeros(1,symbols_per_carrier*(IFFT_bin_length+GI)+GIP);
%注意并串转换后数据的长度为symbols_per_carrier*(IFFT_bin_length
%+GI)+GIP,这里考虑了循环前缀与循环后缀的重叠相加
windowed_Tx_data(1:IFFT_bin_length+GI+GIP)=windowed_time_wave_matrix_ cp(1,:);
%符第一个加窗带循环前缀后缀的OFDM符号值windowed_Tx_data,即发送串
行数据
for i=1:symbols_per_carrier-1;
windowed_Tx_data((IFFT_bin_length+GI)*i+1:(IFFT_bin_length+GI)*(i+1)+GIP)=wi ndowed_time_wave_matrix_cp(i+1,:);
%并串转换，循环前缀与循环后缀重叠相加
end
Tx_data_withoutwindow
=reshape(time_wave_matrix_cp',(symbols_per_carrier)*(IFFT_bin_length+GI+GIP),1)' ;
%不加窗数据并串转换
Tx_data=reshape(windowed_time_wave_matrix_cp',(symbols_per_carrier)*(IFFT _bin_length+GI+GIP),1)';
%加窗数据，单按照循环前缀与循环后缀不重叠相加进行并串转换。此时数据长度为(symbols_per_carrier)*(IFFT_bin_length+GI+GIP)
temp_time1=(symbols_per_carrier)*(IFFT_bin_length+GI+GIP);
%加窗，循环前缀与循环后缀不重叠数据长度，纪委发送的总数据比特数
figure(5);
subplot(2,1,1);
plot(0:temp_time1-1,Tx_data);
%画循环前缀与循环后缀不顾从谍相加OFDM信号的时域波形。
grid on
ylabel('Amplitude(volts)')
xlabel('Time(samples)')
title('OFDM Time Signal')
temp_time2=symbols_per_carrier*(IFFT_bin_length+GI)+GIP;
%加窗，循环前缀与循环后缀重叠相加数据长度
subplot(2,1,2);
plot(0:temp_time2-1,windowed_Tx_data);
%画循环前缀与循环后缀重叠相加OFDM信号的时域波形
grid on
ylabel('Amplitude(volts)')
xlabel('Time (samples)')
title('OFDM Time Signal')
%未加窗发送信号频谱
%对发送数据分段，分段计算频谱，取平均值作为OFDM信号的频谱
symbols_per_average=ceil(symbols_per_carrier/5);
avg_temp_time=(IFFT_bin_length+GI+GIP)*symbols_per_average;
averages=floor(temp_time1/avg_temp_time);
%将发送数据分5段，每段数据长度为avg_temp_time
average_fft(1:avg_temp_time)=0;%存放平均后的OFDM信号的谱，先置0 for a=0:(averages-1)
subset_ofdm=Tx_data_withoutwindow(((a*avg_temp_time)+1):((a+1)*avg_temp_time ));%分段
subset_ofdm_f=abs(fft(subset_ofdm));%对分段后的OFDM信号计算频谱
average_fft=average_fft+(subset_ofdm_f/averages);%取平均
end
average_fft_log=20*log10(average_fft);%求对数平均谱
figure(6);
subplot(2,1,1)
plot((0:(avg_temp_time-1))/avg_temp_time,average_fft_log);
%画未加窗OFDM符号对数平均谱
hold on
grid on
axis([0 0.5 -20 max(average_fft_log)])
ylabel('Magnitude(dB)')
xlabel('Normalized Frequency(0.5=fs/2)')
title('OFDM Signal Specturm')
%计算加窗OFDM信号的频谱
%这段程序功能与上段程序类似，不同之处在于这段程序是对加窗OFDM信号进行分段计算频谱，再取平均
symbols_per_average=ceil(symbols_per_carrier/5);
avg_temp_time=(IFFT_bin_length+GI+GIP)*symbols_per_average;
averages=floor(temp_time1/avg_temp_time);
%将发送数据分5段，每段数据长度为avg_temp_time
average_fft(1:avg_temp_time)=0;%存放平均后的OFDM信号的谱，先置0
for a=0:(averages-1)
subset_ofdm=Tx_data(((a*avg_temp_time)+1):((a+1)*avg_temp_time));%分段
subset_ofdm_f=abs(fft(subset_ofdm));%对分段后的OFDM信号计算频谱
average_fft=average_fft+(subset_ofdm_f/averages);%取平均
end
average_fft_log=20*log10(average_fft);%求对数平均谱
subplot(2,1,2)
plot((0:(avg_temp_time-1))/avg_temp_time,average_fft_log);
%画加窗OFDM符号对数平均谱
hold on
grid on
axis([0 0.5 -20 max(average_fft_log)])
ylabel('Magnitude(dB)')
xlabel('Normalized Frequency(0.5=fs/2)')
title('Windowed OFDM Signal Specturm')
%=================经过加性高斯白噪声信道
Tx_signal_power=var(windowed_Tx_data);%计算信号功率
linear_SNR=10^(SNR/10);%转换对数信噪比为线性幅度值
noise_sigma=Tx_signal_power/linear_SNR;%计算噪声功率，也就是方差
noise_scale_factor=sqrt(noise_sigma);%计算标准差
noise=randn(1,((symbols_per_carrier)*(IFFT_bin_length+GI))+GIP)*noise_scale_ factor;%产生功率为noise_scale_factor高斯噪声
Rx_data=windowed_Tx_data+noise;
%在发送数据上加噪声，相当于OFDM信号经过加性高斯白噪声信道
%==================OFDM解调======================
Rx_data_matrix=zeros(symbols_per_carrier,IFFT_bin_length+GI+GIP);
%存放并行的接受数据
for i=1:symbols_per_carrier;
Rx_data_matrix(i,:)=Rx_data(1,(i-1)*(IFFT_bin_length+GI)+1:i*(IFFT_bin_length+GI )+GIP);%串并转换
end
Rx_data_complex_matrix=Rx_data_matrix(:,GI+1:GI+IFFT_bin_length);
%去掉循环前缀与循环后缀，取出OFDM符号传输的数据
Y1=fft(Rx_data_complex_matrix,IFFT_bin_length,2);
%求取FFT即OFDM信号解调
Rx_carriers=Y1(:,carriers);
%取出carriers序号对应的子载波上的发送数据，去掉加入的零以及共轭对称部分
Rx_phase=angle(Rx_carriers);
%计算接收信号的相位特性
Rx_mag=abs(Rx_carriers);
%计算接收信号的幅度
[M,N]=pol2cart(Rx_phase,Rx_mag);%转换极坐标数据为直角坐标数据
Rx_complex_carrier_matrix=complex(M,N);%两个直角坐标的实数据为构成复数句
figure(7);
plot(Rx_complex_carrier_matrix,'*r');%画接收信号的星座图
axis([-4,4,-4,4]);
title('SNR=30dB接收信号星座图');
grid on
%16QAM解调
Rx_serial_complex_symbols=reshape(Rx_complex_carrier_matrix',size(Rx_compl ex_carrier_matrix,1)*size(Rx_complex_carrier_matrix,2),1)';
%将矩阵Rx_complex_carrier_matrix转换为1的数组
Rx_decoded_binary_symbols=demoduqam16(Rx_serial_complex_symbols);
%进行16QAM解调
baseband_in=Rx_decoded_binary_symbols;
%将解调恢复的二进制信号存放在baseband _in
%误码率计算
bit_errors=find(baseband_in ~= baseband_out);
%将解调恢复的二进制信号与发送二进制信号比较，寻找误码
bit_error_count=size(bit_errors,2);%计算误码个数
ber=bit_error_count/baseband_out_length;%计算误码率