%系統帶寬100MHz=1/Ts=N_sc*子載波間隔
%輸入信号和多徑衰落信道的采樣間隔必須相同。
%(有效不加CP)OFDM符号長度*子載波間隔=1
%每個OFDM符号時長為35.7μs，每時隙14個符号，0.5ms*2時隙*14符号*35.7μs*10子帕=10ms
W=100000000;%系統帶寬Hz
N_sc_Inter=30000;%子載波間隔Hz
N_sc=3276; %系統子載波數（不包括直流載波）number of subcarrier=273*12
N_fft=4096; % FFT 長度【必須2的幂次方】
N_cp=293; % 循環字首長度、Cyclic prefix為1/4*FFT_bin_length 【長CP=352,短CP=288】
M=256; %256QAM調制,調制階數為log2(M)=8(每個符号比特數)，
EBNO=0:1:40; %仿真信噪比Es/N0 = Eb/N0 + 10log10(k),k是每個符号的有用bit數;Es/N0 = snr + 10log10(m),m是符号周期和采樣之間的比率.
%snr = EB/NO +10*log10(rb/w);對于帶通調制信号，例如幅移鍵控ASK、頻移鍵控 PSK 和正交幅度調制 QAM，需要的傳輸帶寬是相應基帶信号的 2 倍，那麼理論上沒有碼間串擾的最大頻譜效率變為 1碼元/s/Hz。
Nd=280; % 每幀包含的OFDM符号數
pilot_Inter=4; %導頻間隔
L=7; %卷積碼限制長度
tblen=6*L; %Viterbi譯碼器回溯深度
Fs=122880000; %4096*30KHz
bits_Per_Symbol=4;%發送每符号含比特數 log2(M)*編碼效率1/2=4
frame=20;%仿真帕數we will simulate 20帕 to generate 'enough' errors 保證足夠的可靠度
rs=1/(33.33*1e-6);
rb=rs*bits_Per_Symbol;%rb=N_sc*bits_Per_Symbol/(35.7*1e-6);
%Eb=sum(P_data.^2)/(N_sc*Nd*bits_Per_Symbol)%信号的總能量除以信号中包含的比特數
%N0=Eb/EBNO
% 基帶資料資料産生
P_data=randi([0 1],1,N_sc*Nd*bits_Per_Symbol);
% 信道編碼（卷積碼、或交織器）
%卷積碼：前向糾錯非線性碼
%交織：使突發錯誤最大限度的分散化
trellis = poly2trellis(7,[133 171]); %(2,1,7)卷積編碼
code_data=convenc(P_data,trellis);%輸入1位會輸出2位。[1,8*N_sc*Nd]
%---------------246QAM調制---------------------------------------
data_temp1= reshape(code_data,log2(M),[])'; %以每組8比特進行分組，M=256，友善進行256QAM調制[N_sc*Nd,8]
data_temp2= bi2de(data_temp1); %二進制轉化為十進制 [N_sc*Nd,1]，友善一個十進制（8個二進制序列）映射為一個複數
modu_data=qammod(data_temp2,M);%256QAM調制[N_sc*Nd,1]
%[15*d,13*d,11*d,9*d,7*d,5*d,3*d,d,-d,-3*d,-5*d,-7*d,-9*d,-11*d,-13*d,-15*d];
%---------------标簽制作---------------------------------------
%label=[real(modu_data(1:(N_sc*Nd/4)));imag(modu_data(1:(N_sc*Nd/4)))].'; %real複數的實部數值,imag求複數的虛部部分
% label=modu_data;
% dlmwrite('256QAM_train_label.csv', label);
%---------------導頻格式---------------------------------------
%pilot_len=Nd/4;%導頻長度=70個十進制調制複數
pilot_len=Nd;%bits_Per_Symbol=4;%每符号含比特數
%pilot_len=Nd*2*frame;%frame=20;%仿真帕數
pilot_symbols=round(rand(1,pilot_len));%導頻符号個數。rand産生1行p列0-1數。round四舍五入取整。=randi([0 1],1,pilot_len);
for i=1:pilot_len
pilot_symbols(1,i)=3+9*1i;
end
%----------------計算導頻和資料數目----------------------------
num_pilot=ceil(N_sc/pilot_Inter);%将3276/4的結果向正無窮方向取整 3276/4=819
if rem(N_sc,pilot_Inter)==0%rem 整除x/y的餘數 2
num_pilot=num_pilot+1;%導頻數=820
end
num_data=N_sc+num_pilot;%資料數=載波個數+導頻數4096
%----------------導頻位置計算----------------------------------
pilot_Indx=zeros(1,num_pilot);%[1,820] 820個導頻
Data_Indx=zeros(1,num_pilot*(pilot_Inter+1));%
for i=1:num_pilot-1
pilot_Indx(1,i)=(i-1)*(pilot_Inter+1)+1;%[1,6,11,16,21,26,...4091]導頻加入的位置 330導頻
end
pilot_Indx(1,num_pilot)=num_data;%為[1,6,11,16,21,26,...4091，4096]
for j=0:num_pilot%0：820
Data_Indx(1,(1+j*pilot_Inter):(j+1)*pilot_Inter)=(2+j*(pilot_Inter+1)):((j+1)*(pilot_Inter+1));
end%第1到4>[2到5] 後延一位。5到8行>[7到10]...3273到3276行>[4092到4095]..3277到3280行>[4097到4100]..3281到3284行>[4102到4105]
Data_Indx=Data_Indx(1,1:N_sc);%取3276個
%Data_Indx與pilot_Indx位置互補為4096
%----------------導頻插入-------------------------------------
piloted_ofdm_syms=zeros(num_data,pilot_len);%[4096,280]
piloted_ofdm_syms(Data_Indx,:)=reshape(modu_data,N_sc,[]);%[3276,280]
piloted_ofdm_syms(pilot_Indx,:)=repmat(pilot_symbols,num_pilot,1);%[820,280]将pilot_symbols整體塊重複為820*1的矩陣
%串并轉換
%modu_data1=reshape(modu_data,N_sc,[]);
%----------------IFFT-------------------------------------
ifft_data=ifft(piloted_ofdm_syms)/sqrt(170); %[4096,280] 歸一化
% 插入保護間隔、循環字首
Tx_cd=[ifft_data(N_fft-N_cp+1:N_fft,:);ifft_data];%把矩陣ifft的末尾N_cp個數補充到最前面%[4389,280] [a;b]縱向拼接上a下b，[a,b]橫向拼接，左右。
% 并串轉換
Tx_data=reshape(Tx_cd,[],1);%由于傳輸需要[4389*280,1]
%-----------% 信道（通過多經瑞利信道、或信号經過AWGN信道）----------
Ber=zeros(1,length(EBNO));
Ser=zeros(1,length(EBNO));
berTheory=zeros(1,length(EBNO));
for jj=1:length(EBNO)
%snr = EBNO(jj)+10*log10(rb/W);%對于“awgn”函數，我們應該使用SNR參數
snr = EBNO(jj)+10*log10(1/2)+10*log10(8);%SNR = EbN0 + 10log10(nBits*coderate) - 10log10(0.5(實數)or1（複數） * upfactor上采樣倍數？)
%單徑信道
%rx_channel=awgn(Tx_data,snr,'measured');%添加高斯白噪聲
%多徑瑞利信道
%chan=rayleighchan(ts,fd,tau,pqb);
%y = filter(chan,Tx_data)；
%fd信道的最大多普勒頻移 機關hz.與速率的換算關系為v×fc/c，fc是載頻 %ts為輸入信号的采樣周期 抽樣間隔機關s.采樣時間=1/fs采樣頻率 %ts=1/fs;%采樣間隔
%tau為每徑相對時延向量 %pdb為每徑相對增益.包含了各徑的功率（當然是均值啦，實際産生的能量都是以此為均值的随機量），以dB為機關
%EPA
%fd=5;
%tau=[0,30,70,90,110,190,410].*1e-9;
%pdb=[0,-1,-2,-3,-8,-17.2,-20.8];
%rx_channel=rayleighchan(ts,fd,tau,pqb);
%y = filter(chan,Tx_data);
%EVA
%fd=70;
%tau=[0,30,150,310,370,710,1090,1730,2510].*1e-9;
%pdb=[0,-1.5,-1.4,-3.6,-0.6,-9.1,-7.0,-12.0,-16.9];
%ETU
%fd=70;
%tau=[0,50,120,200,230,500,1600,2300,5000].*1e-9;
%pdb=[-1.0,-1.0,-1.0,-0,-0,-0,-3.0,-5.0,-7.0];
% SUI
% Ts=1/Fs; %采樣間隔
% Fd=0.5; %Doppler頻偏，Hz
% tau=[0, 0.4e-6, 0.9e-6];
% pdb=[0, -5, -10];
% chan=rayleighchan(Ts, Fd, tau, pdb);%N×M的矩陣，在這裡N是信道處理的資料長（也就是輸入信号的長度），M是多徑數
% channel=filter(chan,Tx_data);%将信道的影響加在輸入的資料上。輸入信号（以Ts為間隔的采樣點）與濾波器進行卷積輸出
% rx_channel=awgn(channel,snr,'measured');%添加高斯白噪聲
Ts=1/Fs; %采樣間隔
Fd=5; %Doppler頻偏，Hz
tau=[0,30,70,90,110,190,410].*1e-9; %多徑延時向量，s
pdb=[0,-1,-2,-3,-8,-17.2,-20.8]; %多徑信道增益向量，dB
chan=rayleighchan(Ts, Fd, tau, pdb);%N×M的矩陣，在這裡N是信道處理的資料長（也就是輸入信号的長度），M是多徑數
channel=filter(chan,Tx_data);%将信道的影響加在輸入的資料上。輸入信号（以Ts為間隔的采樣點）與濾波器進行卷積輸出
rx_channel=awgn(channel,snr,'measured');%添加高斯白噪聲
%%在fading channel下，信道估計的意義就大的多了。LS增加noise的效果
% 串并轉換
Rx_data1=reshape(rx_channel,N_fft+N_cp,[]);%[4389,280]
%---------------% 去掉保護間隔、循環字首---------------------
Rx_data2=Rx_data1(N_cp+1:N_fft+N_cp,:);%[4096,280]
%----------------% FFT-------------------------------------
Rx_carriers=sqrt(170)*fft(Rx_data2);%[4096,280] 歸一化
%----------------導頻和資料提取--------------------------------
Rx_pilot=Rx_carriers(pilot_Indx,:);%導頻提取
Rx_fre_data=Rx_carriers(Data_Indx,:);%資料提取
%----------------% 信道估計與插值（均衡)---------------------------------
pilot_patt=repmat(pilot_symbols,num_pilot,1);
h=Rx_pilot./pilot_patt;
H=interp1( pilot_Indx(1:end),h,Data_Indx(1:end),'linear','extrap');%分段線性插值：插值點處函數值由連接配接其最鄰近的兩側點的線性函數預測。對超出已知點集的插值點用指定插值方法計算函數值
% 信道校正
output=Rx_carriers(Data_Indx(1:end),:)./H;%[3300,70]
%---------------星座圖---------------------------------------
% output_ls=reshape(output,1,[]);
% figure(2);
% scatterplot(output_ls),grid;
%--------------% 資料制作---------------------------------------
data_IQ =[real(Rx_fre_data(1:(N_sc*Nd)));imag(Rx_fre_data(1:(N_sc*Nd)))].'; %real複數的實部數值,imag求複數的虛部部分
label = zeros(N_sc*Nd, M, 'int8');
for i= 1: (N_sc*Nd)
label(i,data_temp2(i)+1)=1; %256種
end
dlmwrite('train_Data256QAM.csv', data_IQ);
dlmwrite('train_label256QAM.csv', label);
%RNN訓練後csv資料讀取
% filename='after_train_data.csv'；
%M = csvread(filename) 輸入檔案名，将資料讀到矩陣M中
%De_data1 = reshape(M,[],1);
%----------------% 256QAM解調------------------------------------
demodulation_data=qamdemod(output,M);
De_data1 = reshape(demodulation_data,[],1);
De_data2 = de2bi(De_data1);%十進制轉換為二進制
De_Bit = reshape(De_data2',1,[]);
%----------------% （解交織）-------------------------------------
% 信道譯碼（維特比譯碼）
trellis = poly2trellis(7,[133 171]);
rx_c_de = vitdec(De_Bit,trellis,tblen,'trunc','hard'); %硬判決
%----------------% 計算誤碼率------------------------------------
%[err,Ber2(jj)] = biterr(De_Bit(1:length(code_data)),code_data);%譯碼前的誤碼率[number錯誤比特的個數,ratio誤比特率] = biterr(xin,yout)
%[err, Ber(jj)] = biterr(rx_c_de(1:length(P_data)),P_data);%譯碼後的誤碼率
%[err,Ber2(jj)] = biterr(code_data,De_Bit(1:length(code_data)));%譯碼前的誤碼率[number錯誤比特的個數,ratio誤比特率] = biterr(xin,yout)
[number_err,Ber(jj)] = biterr(P_data,rx_c_de(1:length(P_data)));%譯碼後的誤碼率
% error_symbol_data1= reshape(rx_c_de,log2(M),[])'; %以每組8比特進行分組，輸出兩列資料
% error_symbol_data_receive=bi2de(error_symbol_data1);
% error_symbol_data2= reshape(P_data,log2(M),[])'; %以每組8比特進行分組，輸出兩列資料
% error_symbol_data_transmite=bi2de(error_symbol_data2);%輸出一列資料
% [err,Ser(jj)] = symerr(error_symbol_data_transmite,error_symbol_data_receive);
berTheory = berawgn(EBNO,'qam',M);%ber = berawgn(EbNo,'qam',M)EbNo是比特能量與噪聲功率頻譜密度的比率，以dB為機關。
%作為該berawgn功能的替代方法，調用BERTool GUI（bertool），然後使用“ 理論”頁籤。
%Bit error rate (BER) for uncoded AWGN channels
%trellis = poly2trellis(7,[133 171]);
%spect = distspec(trellis,2)
%bercoding(EBNO,'conv','hard',1/2,spect)%coded AWGN channels
end
%berfit()
figure(3);
%semilogy(EBN0,BER0,'b-s');
%plot(Rayebn0,Rayber0,'b-s');
%plot(SNR,Ser,'b-s');x
%hold on;
%semilogy(SNR,Ber,'r-o');
plot(EBNO,Ber,'r-o');
hold on;
legend('256QAM調制、基于導頻LS');
%legend('bertool256QAM調制、rayleigh、分級階數7','LS估計256QAM調制、rayleigh多徑+AWGN、卷積碼');
hold on;
xlabel('Eb/N0');
ylabel('BER');
axis([0 40 0 1]);
grid on;
title('誤比特率曲線');
%title('AWGN下誤比特率曲線');
% figure(4)
% subplot(2,1,1);
% x=0:1:50;
% stem(x,P_data(1:51));
% ylabel('振幅');
% title('發送資料（前50個資料為例)');
% legend('256QAM調制、卷積譯碼');
% subplot(2,1,2);
% x=0:1:50;
% stem(x,rx_c_de(1:51));
% ylabel('振幅');
% title('接收資料(以前50個資料為例)');
% legend('256QAM調制、卷積譯碼');
%a=3*log2(L);b=L^(2)-1;c=1-L^(-1);
%theo_sym_rate(n)=(erfc(sqrt((10.^(snrpBit(n).*0.1)).*a./b)).*2*c).*(1-erfc(sqrt((10.^(snrpBit(n).*0.1)).*a./b)).*(1/2)*c);%計算理論誤碼率
![數學模組化淺談[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)