一种使用16QAM的OFDM系统MATLAB仿真

最新推荐文章于 2024-06-27 08:11:34 发布

DanielWu2020

最新推荐文章于 2024-06-27 08:11:34 发布

阅读量3.5k

点赞数 14

文章标签： matlab 经验分享其他

本文链接：https://blog.csdn.net/DanielWu2020/article/details/121717197

版权

分享一种完整的OFDM系统的MATLAB程序，供大家学习参考。

本程序采用16QAM编码，并详细解释了各个参数的值。

本代码总共分为四个函数部分：主函数、16QAM调制函数、解调函数以及窗函数模块。

1、主函数：

%传输速率 25Mbit/s
%可容忍的时延扩展 200ns
%带宽<18MHz
%200ns的时延扩展意味着保护间隔的有效取值为800ns(一般为时延扩展均方根值的2~4倍)
%OFDM符号周期长度为保护间隔的6倍，6X800ns=4.8us
%子载波间隔取4.8us-0.8us=4us的倒数，即250kHz
%每个符号需要传送(25Mbit/s)/(1/4.8us) = 120bit
%若用16QAM和码率为1/2的编码方法，每个子载波可以携带2bit有效信息，需要用60个子载波来满足每个符号120bit的传输速率
%希望在FFT/IFFT间隔内正好有64个采样点，以保持子载波之间的正交性。采样速率为64/4us=16MHz
%然而在4.8us内，这一采样速率不能保证采样数量为整数，因此将每符号的采样速率设定为78，即采样速率为78/4.8=16.25MHz
%FFT/IFFT运算时间产犊为64/16.25MHz=3.9385us
%此时保护间隔和子载波间隔都稍微大于FFT/IFFT运算时间长度为4us时的情况，分别为861.5us和253kHz
clear all;
close all;
carrier_count = 200;
%OFDM子载波数为512，其中400（2*carrier_count）为数据符号，其余赋0值
symbols_per_carrier = 20;
%每个子载波上的符号数，这里即为OFDM符号的个数
bits_per_symbol = 4;
%OFDM符号的每个子载波上传输的比特数，4bit通常采用16QAM
IFFT_bin_length = 512;%FFT长度，一个OFDM符号的子载波个数为512
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 = 30;%信噪比30dB
%----------------OFDM信号产生-------------------
baseband_out_length = carrier_count * symbols_per_carrier * bits_per_symbol;
%计算传输数据总的比特数，为200x20x4bit=16000bit。16000bit的构成为20个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;
%计算OFDM符号子载波的序号，conjugate_carriers中存放的序号是282~481
rand('twister',0);
baseband_out = round(rand(1,baseband_out_length));
%baseband_out产生16000bit待传输的二进制比特流，存放发送的二进制信号与后面借条后的二进制信号比较，剋以计算误码率。
complex_carrier_matrix = qam16(baseband_out);
%16QAM调制并绘制星座图
complex_carrier_matrix = reshape...
    (complex_carrier_matrix',carrier_count,symbols_per_carrier)';
%转换complex_carrier_matrix中的数据为carrier_count,symbols_per_carrier矩阵，20x200矩阵
figure(1);
plot(complex_carrier_matrix,'*r');%
axis([-4,4,-4,4]);
title('16QAM调制后星座图');
grid on;
%IFFT，OFDM调制
IFFT_modulation = zeros(symbols_per_carrier,IFFT_bin_length);
%将symbols_per_carrier*IFFT_bin_length的矩阵赋0，20x512矩阵，512是IFFT长度，也是OFDM符号的子载波个数
IFFT_modulation(:,carriers) = complex_carrier_matrix;%XX(:,1:5)
%将20x200的complex_carrier_matrix数据赋给IFFT_modulation的第29~228列
IFFT_modulation(:,conjugate_carriers) = conj(complex_carrier_matrix);
%将20x200的complex_carrier_matrix数据赋给512个子载波的第282~481个子载波，这样仅有400个子载波有数据，其余为0，相当于补0
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,:));
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('Amplitude');
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(beta,IFFT_bin_length + GI)';
    %调用rcoswindow产生升余弦窗，加窗
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]);
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) = windowed_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);
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);
grid on;
ylabel('Amplitude(volts)');
xlabel('Time(samples)');
title('OFDM Time Signal');
%画循环前缀循环后缀重叠相加的时域波形
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
averages_fft(1:avg_temp_time) = 0;%存放平均后的OFDM信号的谱，先置零
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信号计算频谱
    averages_fft = averages_fft + (subset_ofdm_f/averages);%取平均
end
averages_fft_log = 20 * log10(averages_fft);%求对数平均谱
figure(6);
subplot(2,1,1);
plot((0:(avg_temp_time - 1))/avg_temp_time,averages_fft_log);
%画未加窗OFDM符号对数平均谱
hold on
grid on
axis([0 0.5 -20 max(averages_fft_log)])
ylabel('Magnitude(dB)');
xlabel('Normalized Frequency(0.5 = fs/2)');
title('OFDM Signal Spectrum');
%计算加窗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);
averages_fft(1:avg_temp_time) = 0;%将数据分成5段，每段数据长度为avg_temp_time
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信号计算频谱
    averages_fft = averages_fft + (subset_ofdm_f/averages);%取平均
end
averages_fft_log = 20 * log10(averages_fft);%求对数平均谱
subplot(2,1,2);
plot((0:(avg_temp_time - 1))/avg_temp_time,averages_fft_log);
%画加窗OFDM符号对数平均谱
hold on
grid on
axis([0 0.5 -20 max(averages_fft_log)])
ylabel('Magnitude(dB)');
xlabel('Normalized Frequency(0.5 = fs/2)');
title('Windowed OFDM Signal Spectrum');
%-----------------------------经过加性噪声信道------------------------
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解调——————————————————
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_data_complex_matrix = complex(M,N);%两个直角坐标的实数据为构成复数据
figure(7)
plot(Rx_data_complex_matrix,'*r');
axis([-4,4,-4,4]);
title('SNR=30dB接受数据星座图');
grid on
%16QAM解调
Rx_serial_complex_symbols = reshape(Rx_data_complex_matrix',...
    size(Rx_data_complex_matrix,1) * size(Rx_data_complex_matrix,2),1)';
%将矩阵Rx_serial_complex_symbols转换为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;%计算误码率

2、16QAM函数模块：

function [complex_qam_data] = qam16(bitdata)
%输入参数：bitdata为二进制数码流
%输出参数：complex_qam_data为16QAM复信号
X1 = reshape(bitdata,4,length(bitdata)/4)';%将二进制数码流以4比特分段
d = 1;
%转换4比特二进制码为十进制码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信号

3、解调模块

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_complex_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_metric decode_symble(i)] = min(metrics);
end
%将接受数据与标准16QAM信号比，找到差最小的，将其对应回复成标准的16QAM信号
decoded_bit_symble = de2bi((decode_symble - 1)','left-msb');
%16QAM转为二进制
demodu_bit_symble = reshape(decoded_bit_symble',...
    1,length(Rx_serial_complex_symbols) * 4);%转换为一行

4、窗函数模块

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:(1 + beta) * Ts + 1
    rcosw(j - 1) = 0.5 + 0.5 * cos((t(j) - Ts) * pi/(beta * Ts));
    %计算升余弦窗第三部分
end
rcosw = rcosw';%转化为列矢量