基於16QAM調製的OFDM系統仿真實現（MATLAB）

時間  2021-01-02
標籤 web 算法 typescript svg spa .net code xml blog 欄目 MATLAB 简体版
原文   原文鏈接

    本仿真系統主要針對光OFDM系統，但須要注意的是，本系統還不可以直接用於採用強度調製直接檢測（IM/DD）的光通訊系統。由於採用IM/DD的光通訊系統，其發送端是把要傳送的信息轉變爲電流信號注入LD或LED，從而得到相應的光信號，這就要求要傳送的信號必須爲正且實的信號，然而經本系統中IFFT獲得的時域信號雖然全爲實的，但含有負數，這對於採用IM/DD的光通訊系統是不行的。
    對此有兩種解決方案，一是在獲得的時域信號上加直流偏置（對應DCO-OFDM），使信號全變爲正值，這樣只須要在接收端去直流偏置便可獲得原信號，可是這種方法浪費發射功率，並且會致使峯值過大，可能使LD或LED工做在非線性區域；二是修改埃爾米特對稱的結構，使全部的偶數位子載波均不傳輸數據（對應ACO-OFDM），這樣經IFFT變換後的數據遵循半波對稱原則，能夠在不丟失任何信息的狀況下去掉複數信號，而後在接收端在相應位置恢復便可，但這種方法雖然克服了第一種方法的缺點，但有一半的子載波不傳輸數據，所以具備較低的頻譜利用率。在實際應用中，將根據實際經權衡後做出合理的選擇。
仿真代碼(MATLAB)以下：web

clc
clear all;
close all;
Nt_carr=256;      %子載波數=FFT點數---256
Np_carr=Nt_carr/2-1; %實際子載波數---127
Sig_per_carr=500; %每子載波含符號數---500
bits_per_symbol=4;      %每符號含比特數,16QAM調製---4
CP_Ratio=1/8;
CP_length=CP_Ratio*Nt_carr;%循環前綴長度
TS_length=Nt_carr+CP_length;%訓練序列長度
SNR=20;

%==================================================
%================信號產生===================================
baseband_out_length=Np_carr * Sig_per_carr * bits_per_symbol;  %所輸入的比特數目  500*256*4
% rand( 'twister',0);
baseband_out=round(rand(1,baseband_out_length));%輸出待調製的二進制比特流

%==============16QAM調製====================================
complex_carrier_matrix=qam16(baseband_out);%列向量

figure;
plot(complex_carrier_matrix,'*r');%16QAM調製後星座圖
title('16QAM調製星座圖');
axis([-5,5,-5,5]);
grid on   %顯示網格線
axis square

%==============串並變換====================================
%矩陣轉置時附加有取共軛
complex_carrier_matrix1=reshape(complex_carrier_matrix',Np_carr,Sig_per_carr)';%串並變換Sig_per_carr*Nt_carr 矩陣
%complex_carrier_matrix1=conj(reshape(complex_carrier_matrix',Nt_carr,Sig_per_carr)')%兩次轉置

%==============埃爾米特映射====================================
carriers=(1:Np_carr)+1;%共軛對稱子載波映射  複數數據對應的IFFT點座標
conjugate_carriers=Nt_carr-carriers+2;%共軛對稱子載波映射  共軛複數對應的IFFT點座標
IFFT_modulation=zeros(Sig_per_carr,Nt_carr);%添0組成IFFT_bin_length IFFT 運算
IFFT_modulation(:,carriers )=complex_carrier_matrix1 ;%未添加導頻信號 ，子載波映射在此處
IFFT_modulation(:,conjugate_carriers )=conj(complex_carrier_matrix1);%共軛複數映射

%=================IFFT===========================
time_wave_matrix=ifft(IFFT_modulation,Nt_carr,2);%OFDM調製 即IFFT行變換
%時域波形矩陣，行爲每載波所含符號數，列IFFT點數，N個子載波映射在其內，每一行即爲一個OFDM符號

%=====================計算峯均比PARR====================================
PAPR=10*log10(Sig_per_carr*Nt_carr*(max(max(time_wave_matrix)))^2/sum(sum(time_wave_matrix.^2)));

%==============降PAPR預操做====================================
b=[exp(pi/4*1i) exp(5*pi/4*1i)];
for i=3:Np_carr
    complex_carrier_matrix1(:,i)=b(1)*complex_carrier_matrix1(:,i);
end
PAPR_temp=zeros(1,4);
for i=1:2
    for j=1:2
        complex_carrier_matrix1_temp=[b(i)*complex_carrier_matrix1(:,1) b(j)*complex_carrier_matrix1(:,2) complex_carrier_matrix1(:,3:Np_carr)];
        carriers=(1:Np_carr)+1;
        conjugate_carriers=Nt_carr-carriers+2;
        IFFT_modulation=zeros(Sig_per_carr,Nt_carr);
        IFFT_modulation(:,carriers )=complex_carrier_matrix1_temp ;
        IFFT_modulation(:,conjugate_carriers )=conj(complex_carrier_matrix1_temp);
        time_wave_matrix=ifft(IFFT_modulation,Nt_carr,2);
        PAPR_temp((i-1)*2+j)=10*log10(Sig_per_carr*Nt_carr*(max(max(time_wave_matrix)))^2/sum(sum(time_wave_matrix.^2)));
    end
end
[~,index]=min(PAPR_temp);
complex_carrier_matrix1_temp=[b(ceil(index/2))*complex_carrier_matrix1(:,1) b(index-2*(ceil(index/2)-1))*complex_carrier_matrix1(:,2) complex_carrier_matrix1(:,3:Np_carr)];
carriers=(1:Np_carr)+1;
conjugate_carriers=Nt_carr-carriers+2;
IFFT_modulation=zeros(Sig_per_carr,Nt_carr);
IFFT_modulation(:,carriers )=complex_carrier_matrix1_temp ;
IFFT_modulation(:,conjugate_carriers )=conj(complex_carrier_matrix1_temp);
time_wave_matrix=ifft(IFFT_modulation,Nt_carr,2);
        
figure;
stem(0:Nt_carr-1, abs(IFFT_modulation(2,1:Nt_carr)),'b*-')%第一個OFDM符號的頻譜
grid on
axis ([0 Nt_carr -0.5 4.5]);
ylabel('Magnitude');
xlabel('IFFT Bin');
title('第一個OFDM符號各頻率對應的幅度（幅頻響應）');

figure;
stem(0:Nt_carr-1, (180/pi)*angle(IFFT_modulation(2,1:Nt_carr)), 'b*-')
hold on
plot(0:Nt_carr-1, (180/pi)*angle(IFFT_modulation(2,1:Nt_carr)), 'go')
axis ([0 Nt_carr -200 +200])
grid on
ylabel('Phase (degrees)')
xlabel('IFFT Bin')
title('第一個OFDM符號各頻率對應的相位（相頻響應）')

figure;
subplot(211)
plot(0:Nt_carr-1,time_wave_matrix(1,:));%第一個OFDM符號的時域波形
axis([0,Nt_carr-1+CP_length, -0.6, 0.6]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('發送端第一個OFDM符號的時域波形');

%=====================添加循環前綴CP====================================
CP=time_wave_matrix(:,Nt_carr-CP_length+(1:CP_length));
time_wave_matrix_add_CP=[CP,time_wave_matrix];
subplot(212)
plot(0:Nt_carr-1+CP_length,time_wave_matrix_add_CP(1,:));%第一個OFDM符號的時域波形
axis([0,Nt_carr-1+CP_length, -0.6, 0.6]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('發送端第一個OFDM符號的時域波形（加CP）');

%=========================並串轉換======================================
time_wave_sequence=reshape(time_wave_matrix_add_CP',(Nt_carr+CP_length)*Sig_per_carr,1)';

%=====================添加訓練序列TS====================================
% TS_0=zeros(1,TS_length);
% BPSKTable=[-1-8i,1+8i];
% TS_0=2*sqrt(65)*randi([0,1],1,TS_length);
% TS_0=5*sqrt(2)*randi([0,1],1,TS_length);%干擾序列平均功率爲TS序列的五倍
% TS_0=2*randi([0,1],1,TS_length)-1;%干擾序列平均功率與TS序列至關
TS_0=sqrt(2)*randi([0,1],1,TS_length);
BPSKTable=[-1,1];
TS_Schmidl_1=BPSKTable(randi([0,1],1,Nt_carr/2)+1);
TS_Schmidl=[TS_Schmidl_1(Nt_carr/2-CP_length+1:Nt_carr/2) TS_Schmidl_1 TS_Schmidl_1];
TS_Minn_1=BPSKTable(randi([0,1],1,Nt_carr/4)+1);
TS_Minn=[-TS_Minn_1(Nt_carr/4-CP_length+1:Nt_carr/4) TS_Minn_1 TS_Minn_1 -TS_Minn_1 -TS_Minn_1];
TS_Park_1=BPSKTable(randi([0,1],1,Nt_carr/4)+1);
TS_Park_2=zeros(1,Nt_carr/4);
for i=1:Nt_carr/4
    TS_Park_2(i)=conj(TS_Park_1(Nt_carr/4-i+1));
end
TS_Park=[conj(TS_Park_2(Nt_carr/4-CP_length+1:Nt_carr/4)) TS_Park_1 TS_Park_2 conj(TS_Park_1) conj(TS_Park_2)];
TS=[TS_0 TS_Schmidl TS_0 TS_Minn TS_0 TS_Park TS_0];
time_wave_sequence_add_TS=[TS time_wave_sequence];

figure;
subplot(311)
plot(0:Nt_carr-1+CP_length,TS_Schmidl);%訓練序列波形
axis([0,Nt_carr-1+CP_length, -1.5, 1.5]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('Schmidl訓練序列波形');

subplot(312)
plot(0:Nt_carr-1+CP_length,TS_Minn);%訓練序列波形
axis([0,Nt_carr-1+CP_length, -1.5, 1.5]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('Minn訓練序列波形');

subplot(313)
plot(0:Nt_carr-1+CP_length,TS_Park);%訓練序列波形
axis([0,Nt_carr-1+CP_length, -1.5, 1.5]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('Park訓練序列波形');

%=====================通過高斯噪聲信道====================================
figure;
received_time_wave_sequence_add_TS=awgn(time_wave_sequence_add_TS,SNR,'measured');
received_TS=received_time_wave_sequence_add_TS(1:(7*TS_length));
received_TS_Schmidl=received_TS((TS_length+1):(2*TS_length));
received_TS_Minn=received_TS((3*TS_length+1):(4*TS_length));
received_TS_Park=received_TS((5*TS_length+1):(6*TS_length));
subplot(311)
plot(0:Nt_carr-1+CP_length,received_TS_Schmidl);%接收端Schmidl訓練序列波形
axis([0,Nt_carr-1+CP_length, -1.5, 1.5]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('接收端Schmidl訓練序列波形');

subplot(312)
plot(0:Nt_carr-1+CP_length,received_TS_Minn);%接收端Minn訓練序列波形
axis([0,Nt_carr-1+CP_length, -1.5, 1.5]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('接收端Minn訓練序列波形');

subplot(313)
plot(0:Nt_carr-1+CP_length,received_TS_Park);%接收端Park訓練序列波形
axis([0,Nt_carr-1+CP_length, -1.5, 1.5]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('接收端Park訓練序列波形');

%=====================Schmidl方法進行符號定時====================================
Syn_length=3*TS_length/2;
p_Schmidl=zeros(1,Syn_length);
R_Schmidl=zeros(1,Syn_length);
M_Schmidl=zeros(1,Syn_length);
Syn_Schmidl=received_TS(1:(3*TS_length));
for m=1+CP_length:Syn_length+CP_length
    for k=0:Nt_carr/2-1
        p_Schmidl(m-CP_length)=p_Schmidl(m-CP_length)+conj(Syn_Schmidl(m+k))*Syn_Schmidl(m+k+Nt_carr/2);
        R_Schmidl(m-CP_length)=R_Schmidl(m-CP_length)+abs(Syn_Schmidl(m+k+Nt_carr/2))^2;
    end
    M_Schmidl(m-CP_length)=abs(p_Schmidl(m-CP_length))^2/R_Schmidl(m-CP_length)^2;
end
[maximum_Schmidl,max_index_Schmidl]=max(M_Schmidl);

%=====================Minn方法進行符號定時====================================
p_Minn=zeros(1,Syn_length);
R_Minn=zeros(1,Syn_length);
M_Minn=zeros(1,Syn_length);
Syn_Minn=received_TS((2*TS_length+1):(5*TS_length));
for m=CP_length+1:Syn_length+CP_length
    for l=0:1
        for k=0:Nt_carr/4-1
            p_Minn(m-CP_length)=p_Minn(m-CP_length)+conj(Syn_Minn(m+k+l*Nt_carr/2))*Syn_Minn(m+k+l*Nt_carr/2+Nt_carr/4);
            R_Minn(m-CP_length)=R_Minn(m-CP_length)+abs(Syn_Minn(m+k+l*Nt_carr/2+Nt_carr/4))^2;
        end
    end
    M_Minn(m-CP_length)=abs(p_Minn(m-CP_length))^2/R_Minn(m-CP_length)^2;
end
[maximum_Minn, max_index_Minn]=max(M_Minn);

%=====================Park方法進行符號定時====================================
p_Park=zeros(1,Syn_length);
R_Park=zeros(1,Syn_length);
M_Park=zeros(1,Syn_length);
Syn_Park=received_TS((4*TS_length+1):(7*TS_length));
for m=CP_length+1:3*TS_length/2+CP_length
    for k=0:Nt_carr/2-1
        p_Park(m-CP_length)=p_Park(m-CP_length)+Syn_Park(m+k)*Syn_Park(m-k+Nt_carr-1);
        R_Park(m-CP_length)=R_Park(m-CP_length)+abs(Syn_Park(m+k))^2;
    end
    M_Park(m-CP_length)=abs(p_Park(m-CP_length))^2/R_Park(m-CP_length)^2;
end
[maximum_Park, max_index_Park]=max(M_Park);

figure;
plot(1:Syn_length,M_Schmidl,'r')
hold on
plot(1:Syn_length,M_Minn,'g')
hold on
plot(1:Syn_length,M_Park,'b')
grid on; 
% axis([0,400,0,1.1]); 
title('三種方法獲得的定時量度曲線'); 
legend('Schmidl算法','Minn算法','Park算法')
xlabel('Time(sample)'); 
ylabel('Timing Metric');

%===========================串並轉換========================================
received_time_wave_matrix_add_CP=reshape(received_time_wave_sequence_add_TS((7*TS_length+1):end)',Nt_carr+CP_length,Sig_per_carr)';%去除訓練序列
figure;
subplot(311)
plot(0:Nt_carr-1+CP_length,received_time_wave_matrix_add_CP(1,:));%接收端第一個OFDM符號的時域波形（含CP）
axis([0,Nt_carr-1+CP_length, -0.6, 0.6]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('接收端第一個OFDM符號的時域波形（含CP）');

%=====================去除循環前綴CP====================================
received_time_wave_matrix=received_time_wave_matrix_add_CP(:,CP_length+1:CP_length+Nt_carr);%去CP
subplot(312)
plot(0:Nt_carr-1,real(received_time_wave_matrix(1,:)));%接收端第一個OFDM符號的時域波形（含CP）
axis([0,Nt_carr-1+CP_length, -0.6, 0.6]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('接收端去除CP後的第一個OFDM符號的時域波形');

%=====================信道估計與均衡====================================
% TS_FFT=fft(TS_Schmidl((CP_length+1):TS_length));
% received_TS_FFT=fft(received_TS_Schmidl((CP_length+1):TS_length));
% H=received_TS_FFT./TS_FFT;

H=received_TS_Schmidl((CP_length+1):TS_length)./TS_Schmidl((CP_length+1):TS_length);

% H=fft(received_time_wave_matrix(1,:))./fft(time_wave_matrix(1,:));

received_time_wave_matrix_FFT=fft(received_time_wave_matrix,Nt_carr,2);
received_time_wave_matrix_equilibrium_FFT=received_time_wave_matrix_FFT./H;

received_time_wave_matrix_equilibrium=ifft(received_time_wave_matrix_equilibrium_FFT,Nt_carr,2);
subplot(313)
plot(0:Nt_carr-1,real(received_time_wave_matrix_equilibrium(1,:)));%接收端第一個OFDM符號的時域波形（含CP）
axis([0,Nt_carr-1+CP_length, -0.6, 0.6]);
grid on;
ylabel('Amplitude');
xlabel('Time');
title('接收端均衡後第一個OFDM符號的時域波形');

%=================================FFT======================================
received_equilibrium_Hermite=fft(received_time_wave_matrix_equilibrium,Nt_carr,2);

%===========================解埃爾米特映射==================================
received_equilibrium=received_equilibrium_Hermite(:,carriers);

%===========================下降PAPR後處理==================================
received_equilibrium(:,3:Np_carr)=received_equilibrium(:,3:Np_carr)/b(1);
received_equilibrium(:,1)=received_equilibrium(:,1)/b(ceil(index/2));
received_equilibrium(:,2)=received_equilibrium(:,2)/b(index-2*(ceil(index/2)-1));

%=============================並串轉換=====================================
received_complex_carrier_matrix1=reshape(received_equilibrium',Np_carr*Sig_per_carr,1)';

figure;
plot(received_complex_carrier_matrix1,'*r');%接收端星座圖
title('接收端(16QAM解調製前)星座圖');
axis([-5,5,-5,5]);
grid on   %顯示網格線
axis square

%===========================16QAM解調==================================
demodu_baseband_out=deqam16(received_complex_carrier_matrix1);
[~,ber]=symerr(demodu_baseband_out,baseband_out);

ber_carriers=zeros(1,Np_carr);
for j=1:Np_carr
    for i=1:Sig_per_carr
        for k=1:bits_per_symbol
            if demodu_baseband_out((i-1)*Np_carr*bits_per_symbol+(j-1)*bits_per_symbol+k)~=baseband_out((i-1)*Np_carr*bits_per_symbol+(j-1)*bits_per_symbol+k)
                ber_carriers(j)=ber_carriers(j)+1;
            end
        end
    end
end
figure;
plot(1:Np_carr,ber_carriers,'--r*');
title('各子載波誤碼個數')
ylabel('誤碼個數');
xlabel('子載波編號');
xlim([1 Np_carr]);
grid on;

OFDM框圖：
在這裏插入圖片描述算法

OFDM信號：
在這裏插入圖片描述在這裏插入圖片描述
在這裏插入圖片描述
星座圖：
在這裏插入圖片描述typescript

請添加圖片描述
請添加圖片描述
同步效果：
請添加圖片描述
請添加圖片描述
更多內容（誤碼率曲線、降PAPR等）請關注：https://blog.csdn.net/hyl1181/article/details/107490533
所有MATLAB文件詳見：https://download.csdn.net/download/hyl1181/12646924svg

本文同步分享在 博客「使君杭千秋」（CSDN）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。ui

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程