基于MATLAB的IS-95系統(tǒng)仿真

1、CDMA 期 末 論 文基于MATLAB的IS-95系統(tǒng)仿真 課程名稱 CDMA蜂窩移動通信 姓 名 學 號 專 業(yè) 成 績 任課老師 上課學期 1 作業(yè)要求l 分析給出的源程序，逐行給出中文注釋，并分析仿真結(jié)果。l 按照如圖1所示的框圖設計一個CDMA系統(tǒng)，并進行仿真。信源卷積編碼交織編碼擴頻加擾碼調(diào)制AWGN信道解調(diào)解擾碼信宿解卷積解交織解擴圖1 CDMA系統(tǒng)框圖2 題目一2.1 源程序注釋2.1.1 主程序部分%main_CDMA_forward.m%此函數(shù)用于IS-95前向鏈路系統(tǒng)的仿真，包括擴%頻調(diào)制，匹配濾波，RAKE接收等相關(guān)通信模塊。%但沒有在頻帶上進行調(diào)制仿真。%仿真環(huán)境:

2、加性高斯白噪聲信道.%數(shù)據(jù)速率 = 9600 KBps% clear allclose allclc disp('-start-');% 定義通信仿真全局變量global Zi Zq Zs show R Gi Gq% Zi I支路擴頻PN序列生成器狀態(tài)% ZQ Q支路擴頻PN序列生成器狀態(tài)% Zs 擾碼PN序列生成器狀態(tài)% show 指示是否畫出仿真過程各部分產(chǎn)生的輸出% R 仿真時碼元寬度% Gi I支路擴頻PN序列% Gq Q支路擴頻PN序列 clear j;show = 0;SD = 0; % 選擇軟/硬判決接收 %-主要的仿真參數(shù)設置-% 信道數(shù)據(jù)率為9600 bpsB

3、itRate = 9600;% 碼片速率為1.2288 MHzChipRate = 1228800;% 消息碼元個數(shù)N = 184;% 匹配濾波器類型-升余弦MFType = 1; R = 5; %+Viterbi生成多項式+% 維特比譯碼器（譯卷積碼）G_Vit = 1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1;% 生成矩陣的長度K = size(G_Vit, 2);% 輸出碼片數(shù)L = size(G_Vit, 1);%+ %+Walsh矩陣+% 規(guī)定需要的Walsh序列長度WLen = 64;% 生成 64 bit Walsh序列1010.Walsh = res

4、hape(1;0*ones(1, WLen/2), WLen , 1);%Walsh = zeros(WLen ,1);%+ %+擴頻調(diào)制PN碼的生成多項式+%Gi = 1 0 1 0 0 0 1 1 1 0 1 0 0 0 0 1'%Gq = 1 0 0 1 1 1 0 0 0 1 1 1 1 0 0 1' % I信道擴頻PN序列生成多項式Gi_ind = 15, 13, 9, 8, 7, 5, 0'% I信道擴頻PN序列生成多項式Gq_ind = 15, 12, 11, 10, 6, 5, 4, 3, 0'% 初始化序列Gi（16 bit，全零）Gi = z

5、eros(16, 1);% 根據(jù)生成多項式Gi_ind將對應位置1Gi(16-Gi_ind) = ones(size(Gi_ind);% 初始化I路信道PN碼生成器的初始狀態(tài)Zi = zeros(length(Gi)-1, 1); 1;% 初始化序列Gi（16 bit，全零）Gq = zeros(16, 1);% 根據(jù)生成多項式Gq_ind將對應位置1,得到二進制生成多項式Gq(16-Gq_ind) = ones(size(Gq_ind);% Q路信道PN碼生成器的初始狀態(tài)Zq = zeros(length(Gq)-1, 1); 1;%+ %+擾碼生成多項式+% 擾碼序列生成多項式Gs_ind

6、 = 42, 35, 33, 31, 27, 26, 25, 22, 21, 19, 18, 17, 16, 10, 7, 6, 5, 3, 2, 1, 0'% 初始化序列GsGs = zeros(43, 1);% 根據(jù)生成多項式Gs_ind將對應位置1,得到二進制生成多項式Gs(43-Gs_ind) = ones(size(Gs_ind);% 長序列生成器的初始狀態(tài)Zs = zeros(length(Gs)-1, 1); 1;%+ %+AWGN信道+% 擴頻增益EbEc = 10*log10(ChipRate/BitRate);%EbEcVit = 10*log10(L);% 仿真信

7、噪比范圍(dB)% EbNo = -2 : 0.5 : 6.5;EbNo = -2 : 0.5 : 4;%+ %- %-主程序- % 初始化向量ErrorsB = ; ErrorsC = ; NN = ;% ErrorsB% ErrorsC% NN% 輸出信息說明采用的維特比譯碼判決類型if (SD = 1) fprintf('n SOFT Decision Viterbi Decodernn');else fprintf('n HARD Decision Viterbi Decodernn');end % 對信噪比范圍內(nèi)每個數(shù)字進行仿真for i=1:leng

8、th(EbNo) % 輸出當前仿真信噪比 fprintf('nProcessing %1.1f (dB)', EbNo(i); iter = 0; % 仿真計數(shù)器 ErrB = 0; % 出錯的比特數(shù)（編碼前） ErrC = 0; % 出錯的子碼數(shù)（編碼后） while (ErrB <300) & (iter <150) % 當誤碼數(shù)超過或者循環(huán)試驗150次以上才可以統(tǒng)計誤碼率 % 激活命令窗口 drawnow; %+發(fā)射機+ % 生成信息碼元（由高斯白噪聲經(jīng)過零比較得到） TxData = (randn(N, 1)>0); % 輸出已編碼、交織并加擾

9、的消息比特與所用到的擾碼 %（碼率：19.2 kbps） TxChips, Scrambler = PacketBuilder(TxData, G_Vit, Gs); % 下面進行擴頻調(diào)制，調(diào)制后碼率為1.2288Mcps x PN MF = Modulator(TxChips, MFType, Walsh); %+ %+%信道產(chǎn)生高斯白噪聲 noise = 1/sqrt(2)*sqrt(R/2)*( randn(size(x) + j*randn(size(x)*10(-(EbNo(i) - EbEc)/20); %j*randn(size(x)*10(-(EbNo(i) - EbEc)/2

10、0); %向已調(diào)信號中混入加性高斯白噪聲 r = x+noise; %+ %+接收機+ % 軟判決，輸出速率為19.2 Kcps 經(jīng)過卷積編碼的序列 RxSD = Demodulator(r, PN, MF, Walsh); % 對Rake接收機解調(diào)后的信號進行過零比較 RxHD = (RxSD>0); if (SD) % 如果選擇了軟判決則會進行軟判決， % 并進行解擾、解交織維特比譯碼 RxData Metric= ReceiverSD(RxSD, G_Vit, Scrambler); else % 進行硬判決，并進行解擾、解交織維特比譯碼 RxData Metric= Receiv

11、erHD(RxHD, G_Vit, Scrambler); end %+ if(show) % 畫出軟判決輸出結(jié)果 subplot(311); plot(RxSD, '-o'); title('Soft Decisions'); % 將判決結(jié)果與發(fā)送序列異或，顯示出錯碼(維特比譯碼前) subplot(312); plot(xor(TxChips, RxHD), '-o'); title('Chip Errors'); % 將收發(fā)數(shù)據(jù)進行異或，顯示出錯碼 subplot(313); plot(xor(TxData, RxData)

12、, '-o'); title('Data Bit Errors. Metric = ', num2str(Metric); pause; end if(mod(iter, 50)=0) % 如果已循環(huán)次數(shù)為50的倍數(shù)則打印.將變量保存 fprintf('.'); % 將變量保存到硬盤 save TempResults ErrB ErrC N iter end % 計算編碼前，解碼后出錯的消息比特(9600 Bps)數(shù)量 ErrB = ErrB + sum(xor(RxData, TxData); % 計算經(jīng)過編碼出錯的子碼(19.2 kBps)數(shù)

13、 ErrC = ErrC + sum(xor(RxHD, TxChips); % 循環(huán)次數(shù)加1 iter = iter+ 1; end % 生成誤比特數(shù)向量 ErrorsB = ErrorsB; ErrB; % 生成誤碼數(shù)向量 ErrorsC = ErrorsC; ErrC; % 計算本次共仿真的消息比特數(shù) NN = NN; N*iter; % 保存所有產(chǎn)生的數(shù)據(jù) save SimData *end %+誤碼率計算+% 計算誤碼率（卷積編碼前的信號）PerrB = ErrorsB./NN;% PerrB1 = ErrorsB1./NN1;% 計算誤碼率（卷積編碼后的信號）PerrC = Err

14、orsC./NN;% BPSK調(diào)制理論誤碼率（無卷積編碼與擴頻）Pbpsk= 1/2*erfc(sqrt(10.(EbNo/10);% 有維特比譯碼但無擴頻的理論誤碼率PcVit= 1/2*erfc(sqrt(10.(EbNo-EbEcVit)/10);% 有擴頻但無維特比譯碼的理論誤碼率Pc = 1/2*erfc(sqrt(10.(EbNo-EbEc)/10);%+ %+性能仿真顯示+figure;semilogy(EbNo(1:length(PerrB), PerrB, 'b-*');hold on;% %semilogy(EbNo(1:length(PerrB1), Pe

15、rrB1, 'k-o');% hold on;semilogy(EbNo(1:length(PerrC), PerrC, 'b-o');grid on;semilogy(EbNo, Pbpsk, 'b-.');semilogy(EbNo, PcVit, 'k-.x');ylabel('BER');semilogy(EbNo, Pc, 'b-.x');xlabel('信噪比/dB');ylabel('誤碼率');%+disp('-end-');%-2.1

16、.2 PacketBuilder函數(shù)用于信號傳輸之前的卷積編碼和交織編碼和加擾碼% *beginning of file*%PacketBuilder.mfunction ChipsOut, Scrambler = PacketBuilder(DataBits, G, Gs); % %此函數(shù)用于產(chǎn)生IS-95前向鏈路系統(tǒng)的發(fā)送數(shù)據(jù)包 %+variables+% DataBits 發(fā)送數(shù)據(jù)（二進制形式）% G Viterbi編碼生成多項式% Gs 長序列生成多項式（擾碼生成多項式）% ChipsOut 輸入到調(diào)制器的碼序列（二進制形式）% Scrambler 擾碼%+ % 全局變量global

17、Zs % 取得卷積碼的約束長度 K = 9K = size(G, 2); % 取得對于每一位消息碼編碼器的輸出位數(shù) L = 2L = size(G, 1); % 計算擴頻后的子數(shù)%(原始比特：9.6 Kbps-> 卷積編碼后：19.2 Kbps ->擴頻后：1.2288 Mbps) N = 64*L*(length(DataBits)+K-1);% 加入尾比特，并進行卷積編碼chips = VitEnc(G, DataBits; zeros(K-1,1); % 交織編碼% 先將編碼器輸出轉(zhuǎn)為 24行16列矩陣INTERL = reshape(chips, 24, 16); % 再將

18、轉(zhuǎn)置矩陣按串行輸出% 速率=19.2 KBps chips = reshape(INTERL', length(chips), 1); % 產(chǎn)生擾碼 % PNGen 用于產(chǎn)生1.2288 Mc/s 的擾碼LongSeq Zs = PNGen(Gs, Zs, N); % 以下進行分頻，得到19.2 kbps 的擾碼Scrambler = LongSeq(1:64:end); % 對編碼后的比特進行加擾（碼率； 19.2 kbps）ChipsOut = xor(chips, Scrambler); %*end of file*2.1.3 Modulator函數(shù)將信號進行擴頻和調(diào)制之后再發(fā)送

19、%*beginningoffile*%VitEnc.mfunction y = VitEnc(G, x); % 此函數(shù)根據(jù)生成多項式進行Viterbi編碼 % %+variables+% G 生成多項式的矩陣% x 輸入數(shù)據(jù)（二進制形式） % y Viterbi編碼輸出序列%+ % 取得對于每一位消息碼編碼器的輸出位數(shù)K = size(G, 1); % 取得消息比特長度L = length(x); % 進行二維卷積運算 yy = conv2(G, x'); % 去除多余項yy = yy(:, 1:L); % 將輸出矩陣變?yōu)橐恍衴 = reshape(yy,K*L, 1); % 模2運算

20、 y = mod(y, 2); %*endoffile*2.1.4 PNGen函數(shù)根據(jù)PN序列的生成多項式產(chǎn)生PN序列2.1.5 VitEnc函數(shù)進行維比特編碼的函數(shù)%*beginningoffile*%PNGen.mfunction y, Z = PNGen(G, Zin, N);% 此函數(shù)是根據(jù)生成多項式和輸入狀態(tài)產(chǎn)生長度為N的偽隨機序列 %+variables+% G 生成多項式% Zin 移位寄存器初始化% N PN序列長度% y 生成的PN碼序列% Z 移位寄存器的輸出狀態(tài)%+ % 取得生成矩陣的長度L = length(G);% 移位寄存器的初始化Z = Zin;% 初始化輸出序列

21、為全零y = zeros(N, 1);for i=1:N % 輸出最高位 y(i) = Z(L); % 異或 Z = xor(G*Z(L), Z); % 移位 Z = Z(L); Z(1:L-1);end%yy = filter(1, flipud(G), 1; zeros(N-1, 1);%yy = mod(yy, 2);%*end of file*2.1.6 函數(shù)Demodulator對接收到的信號進行解調(diào)和解擴%*beginningoffile*%Modulator.mfunction TxOut, PN, MF = Modulator(chips, MFType, Walsh);%此函

22、數(shù)用于實現(xiàn)IS-95前向鏈路系統(tǒng)的數(shù)據(jù)調(diào)制 %+variables+% chips 發(fā)送的初始數(shù)據(jù)% MFType 成型濾波器的類型選擇% Walsh walsh碼% TxOut 調(diào)制輸出信號序列% PN 用于擴頻調(diào)制的PN碼序列% MF 匹配濾波器參數(shù)%+ % 導入全局變量global Zi Zq show R Gi Gq % 計算輸出字碼數(shù)N = length(chips)*length(Walsh); % 輸入速率 = 19.2 KBps, 輸出速率= 1.2288 Mcps% 先將單極性歸二進制碼變換為雙極性碼% 再與Walsh序列相乘，實現(xiàn)正交調(diào)制tmp = sign(Walsh-1

23、/2)*sign(chips'-1/2);% 串并轉(zhuǎn)換chips = reshape(tmp, prod(size(tmp), 1);% 產(chǎn)生I支路PN擴頻序列（1.2288Mc/s）PNi Zi = PNGen(Gi, Zi, N);% 產(chǎn)生Q支路PN擴頻序列（1.2288Mc/s）PNq Zq = PNGen(Gq, Zq, N);% 將I支路與Q支路在復數(shù)域合并表示PN = sign(PNi-1/2) + j*sign(PNq-1/2);% 同時對同相與正交分量進行擴頻擴頻chips_out = chips.*PN; % 對輸出序列進行0插值（加上信道中傳輸?shù)谋Ｗo間隔）chips

24、 = chips_out, zeros(N, R-1);chips = reshape(chips.' , N*R, 1); %成型濾波器，以下代碼構(gòu)造成型濾波器switch (MFType) case 1 %升余弦濾波器（時域） L = 25; L_2 = floor(L/2); n = -L_2:L_2; B = 0.7; MF = sinc(n/R).*(cos(pi*B*n/R)./(1-(2*B*n/R).2); MF = MF/sqrt(sum(MF.2); case 2 %矩形濾波器 L = R; L_2 = floor(L/2); MF = ones(L, 1); MF

25、 = MF/sqrt(sum(MF.2); case 3 %漢明濾波器 L = R; L_2 = floor(L/2); MF = hamming(L); MF = MF/sqrt(sum(MF.2);end % 轉(zhuǎn)置MF = MF(:); % 利用卷積運算進行基帶濾波仿真TxOut = sqrt(R)*conv(MF, chips)/sqrt(2);% 去除冗余信息TxOut = TxOut(L_2+1: end - L_2); if (show) figure; % 畫出匹配濾波器的時域特性（沖擊響應） subplot(211); plot(MF, '-o'); title

26、('Matched Filter'); % 顯示網(wǎng)格 grid on; subplot(212); % 畫出功率譜 psd(TxOut, 1024, 1e3, 113); title('Spectrum');end%*endoffile*2.1.7 函數(shù)ReceiverHD接收機，硬判決，進行解擾、解卷積和解交織%*beginningoffile*%ReceiverHD.mfunction DataOut, Metric = ReceiverHD(HDchips, G, Scrambler); % % 此函數(shù)用于實現(xiàn)基于Viterbi譯碼的硬判決接收機 %+va

27、riables+% SDchips 硬判決RAKE接收機輸入符號% G Viterbi編碼生成多項式矩陣% Scrambler 擾碼序列% DataOut 接收數(shù)據(jù)（二進制形式） % Metric Viterbi譯碼最佳度量%+ % 如果輸入?yún)?shù)只有一個，則采用默認生成矩陣if (nargin = 1) G = 1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1; end % 速率=19.2 KBps % 解擾HDchips = xor(HDchips, Scrambler); % 解交織INTERL = reshape(HDchips, 16, 24); HDchips

28、 = reshape(INTERL', length(HDchips), 1); % 維特比譯碼DataOut Metric = VitDec(G, HDchips, 1); %*endoffile*2.1.8 函數(shù)VitDec硬判決的維比特譯碼%*beginningoffile*%VitDec.mfunction xx, BestMetric = VitDec(G, y, ZeroTail); % % 此函數(shù)是實現(xiàn)硬判決輸入的Viterbi譯碼 %+variables+% G 生成多項式的矩陣% y 輸入的待譯碼序列 % ZeroTail 判斷是否包含0尾% xx Viterbi譯碼

29、輸出序列 % BestMetric 最后的最佳度量%+ % 輸出碼片數(shù) L = size(G, 1); % 生成多項式長度K= size(G, 2); % 狀態(tài)數(shù)N = 2(K-1); % 最大柵格深度T = length(y)/L; % 生成零矩陣 OutMtrx = zeros(N, 2*L); for s = 1:N % 產(chǎn)生全1的矩陣與十進制轉(zhuǎn)成二進制的矩陣相乘 in0 = ones(L, 1)*0, (dec2bin(s-1), (K-1)-'0'); % 產(chǎn)生全1的矩陣與十進制轉(zhuǎn)成二進制的矩陣相乘 in1 = ones(L, 1)*1, (dec2bin(s-1),

30、 (K-1)-'0'); % 兩矩陣對應位置元素相乘之和對2求余 out0 = mod(sum(G.*in0)'), 2); % 兩矩陣對應位置元素相乘之和對2求余 out1 = mod(sum(G.*in1)'), 2); % 兩個矩陣合并 OutMtrx(s, :) = out0, out1; end % 初始狀態(tài) = 100 PathMet = 0; 100*ones(N-1), 1); % 取出第1列的元素PathMetTemp = PathMet(:,1); % 生成零矩陣Trellis = zeros(N, T); % 往矩陣第1列賦予0、1、2的遞

31、增值Trellis(:,1) = 0 : (N-1)' % 矩陣變形 y = reshape(y, L, length(y)/L); for t = 1:T % 從1到T列逐列取出，并轉(zhuǎn)置 yy = y(:, t)' for s = 0:N/2-1 % 對運算結(jié)果取最小值 B0 ind0 = min( PathMet(1+2*s, 2*s+1) + sum(abs(OutMtrx(1+2*s, 0+1:L) - yy).2); sum(abs(OutMtrx(1+(2*s+1), 0+1:L) - yy).2) ); % 對運算結(jié)果取最小值 B1 ind1 = min( Pat

32、hMet(1+2*s, 2*s+1) + sum(abs(OutMtrx(1+2*s, L+1:L) - yy).2); sum(abs(OutMtrx(1+(2*s+1), L+1:L) - yy).2) ); % 將B0與B1矩陣合并并賦值 PathMetTemp(1+s, s+N/2) = B0; B1; Trellis(1+s, s+N/2, t+1) = 2*s+(ind0-1); 2*s + (ind1-1); end PathMet = PathMetTemp; end %生成零矩陣 xx = zeros(T, 1); %如果包含零尾if (ZeroTail) %最佳標志位置1

33、BestInd = 1; else %否則取矩陣最小值賦值 Mycop, BestInd = min(PathMet); end %最佳度量為矩陣被標志的位 BestMetric = PathMet(BestInd); %將(BestInd-1)/(N/2)結(jié)果向下取整%作為Viterbi譯碼輸出序列最后一位xx(T) = floor(BestInd-1)/(N/2); %取數(shù)組元素賦NextState初值NextState = Trellis(BestInd, (T+1); for t=T:-1:2 %NextState/(N/2)向下取整賦值給xx xx(t-1) = floor(Next

34、State/(N/2); %取數(shù)組元素賦值給NextState NextState = Trellis( (NextState+1), t); end %如果零尾if (ZeroTail) %取其中1到end-K+1位 xx = xx(1:end-K+1); end %*endoffile*2.1.9 函數(shù)ReceiverSD用于軟判決的接收機%*beginningoffile*%ReceiverSD.mfunction DataOut, Metric = ReceiverSD(SDchips, G, Scrambler); % % 此函數(shù)用于實現(xiàn)基于Viterbi譯碼的發(fā)送數(shù)據(jù)的恢復 %+v

35、ariables+% SDchips 軟判決RAKE接收機輸入符號% G Viterbi編碼生成多項式矩陣% Scrambler 擾碼序列% DataOut 接收數(shù)據(jù)（二進制形式） % Metric Viterbi譯碼最佳度量%+ if (nargin = 1) G = 1 1 1 1 0 1 0 1 1; 1 0 1 1 1 0 0 0 1; % 若只給出接收序列參數(shù)，則重新定義生成多項式end % 速率=19.2 KBps SDchips = SDchips.*sign(1/2-Scrambler); % 解擾384bit INTERL = reshape(SDchips, 16, 24)

36、; % 解交織 SDchips = reshape(INTERL', length(SDchips), 1); % 速率=19.2 KBps DataOut Metric = SoftVitDec(G, SDchips, 1); % 卷積解碼，判決%*endoffile* 2.1.10 函數(shù)SoftVitDec用于軟判決的維比特譯碼%*beginningoffile*%SoftVitDec.mfunction xx, BestMetric = SoftVitDec(G, y, ZeroTail); % % 此函數(shù)是實現(xiàn)軟判決輸入的Viterbi譯碼 %+variables+% G 生成

37、多項式的矩陣% y 輸入的待譯碼序列 % ZeroTail 判斷是否包含0尾% xx Viterbi譯碼輸出序列 % BestMetric 最后的最佳度量%+ L = size(G, 1); % 輸出碼片數(shù)K= size(G, 2); % 生成多項式的長度 N = 2(K-1); % 狀態(tài)數(shù) T = length(y)/L; % 最大柵格深度 OutMtrx = zeros(N, 2*L); % 存放256種移位寄存器狀態(tài)所對應的編碼for s = 1:N in0 = ones(L, 1)*0, (dec2bin(s-1), (K-1)-'0'); % 初始狀態(tài)為0時對應移位寄

38、存器狀態(tài)in1 = ones(L, 1)*1, (dec2bin(s-1), (K-1)-'0'); % 初始狀態(tài)為1 out0 = mod(sum(G.*in0)'), 2); % 對應的編碼器輸出 out1 = mod(sum(G.*in1)'), 2); OutMtrx(s, :) = out0, out1; % 存到編碼輸出狀態(tài)end OutMtrx = sign(OutMtrx-1/2); % 轉(zhuǎn)換成雙極性碼 PathMet = 100; zeros(N-1), 1); % 初始狀態(tài) = 100 PathMetTemp = PathMet(:,1);

39、Trellis = zeros(N, T); Trellis(:,1) = 0 : (N-1)' y = reshape(y, L, length(y)/L); % y為2*192矩陣for t = 1:T yy = y(:, t); % 取信碼，每兩個一組for s = 0:N/2-1 % 計算收到序列與所有可能發(fā)送序列的距離 B0 ind0 = max( PathMet(1+2*s, 2*s+1) + OutMtrx(1+2*s, 0+1:L) * yy; OutMtrx(1+(2*s+1), 0+1:L)*yy ); B1 ind1 = max( PathMet(1+2*s, 2*s+1) + OutMtrx(1+2*s, L+1