數(shù)字頻帶系統(tǒng)matlab程序

1、實驗八clear all;close all;echo on;t0=3;%信號時長dt1=0.001 %系統(tǒng)時域采樣間隔fs=1/dt1;%系統(tǒng)采樣頻率df=0.2; %所需的頻率分辨率t=0:dt1:t0;% 產(chǎn)生模擬調(diào)制信號m=ones(1,t0/(3*dt1),-2*ones(1,t0/(3*dt1),zeros(1,t0/(3*dt1)+1);%畫出調(diào)制信號波形及頻譜clffigure(1);subplot(121);plot(t,m(1:length(t);% 畫出調(diào)制信號波形R=4.2;axis(0 t0 -R/2 R/2);xlabel('t');ylabel(&

2、#39;調(diào)制信號');subplot(122);M,m,df1,f=T2F4(m,dt1,df,fs);%求出調(diào)制信號頻譜 plot(f,fftshift(abs(M)% 畫出調(diào)制信號頻譜xlabel('f');ylabel('調(diào)制信號頻譜');axis(-5,5,0,2.5);%A/D%低通抽樣定理t=0:dt1:t0;xt=m;f,xf=T2F(t,xt);dt=0.01;%抽樣信號，抽樣頻率為fsfs=10;sdt=1/fs;st=;for t1=sdt:sdt:t0+sdt st=st m(ceil(t1/dt1);endt1=0:sdt:t0;

3、f1,sf=T2F(t1,st); figure(2);subplot(211);t=0:dt1:t0;plot(t,xt(1:length(t);axis(0 t0 -2.5 2.5);title('原始信號');subplot(212);stem(t1,st);title('抽樣信號');axis(0 t0 -2.5 2.5); %基帶信號形成器figure(3);subplot(211);plot(t,xt(1:length(t);title('原始模擬信號');axis(0 t0 -2.5 2.5); out1=pcm_encode(st

4、);out=2*out1-1;%將編碼結(jié)果由單極性不歸0轉(zhuǎn)換為雙極性不歸0碼dt=0.1;y=;for k=1:length(out) y=y ones(1,1/dt)*out(k);endsubplot(212);t=0:t0/(length(y)-1):t0;plot(t,y);xlabel('t');title('PCM編碼結(jié)果');axis(0 t0 -1.5 1.5); A=1; %載波振幅fc=2; %載波頻率（Hz）snr=15; %信噪比dBN_sample=8;% 基帶信號中每個碼元的的采樣點數(shù)N=length(out1); % 碼元數(shù)為編碼碼

5、元數(shù)Ts=1; % 碼元寬度df=0.01;%頻率分辨率B=1/Ts;f_start=fc-B;f_cutoff=fc+B;fs=fc*N_sample%系統(tǒng)采樣頻率，即考慮載波后，一個碼元內(nèi)的采樣點數(shù)ts=Ts/fs; % 系統(tǒng)采樣間隔t=0:ts:N*Ts-ts;Lt=length(t); d=out;dd=sigexpand(out1,fc*N_sample);gt=ones(1,fc*N_sample); % NRZ波形d_NRZ=conv(dd,gt);d_NRZ1=d_NRZ(1:Lt); d_NRZ1=mabianhuan(d_NRZ1,fs);%進行2DPSK調(diào)制的碼變換，獲得

6、相對碼figure(4);out=2*d_NRZ1-1;%將相對碼由單極性不歸0轉(zhuǎn)換為雙極性不歸0碼y=;dt=0.1; for k=1:length(out) y=y ones(1,1/dt)*out(k);endsubplot(211);t=0:t0/(length(y)-1):t0;plot(t,y);xlabel('t');title('相對碼');axis(0 t0 -1.5 1.5); d_sjx=2*d_NRZ1-1;%生成的雙極性NRZ信號d_sjx1=d_sjx(1:Lt); subplot(212);d_sjx1f,d_sjx1,df1,f=

7、T2F4(d_sjx1,ts,df,fs);%求出雙極性NRZ信號功率譜plot(f,10*log10(abs(fftshift(d_sjx1f).2/length(f);% 畫出雙極性NRZ信號功率譜axis(-3*B 3*B -100 10);xlabel('f');title('相對碼的的雙極性信號PDF'); t=0:ts:N*Ts-ts;%-畫出數(shù)字頻帶信號及其功率譜ht=A*sin(2*pi*fc*t);%載波 % 數(shù)字基帶信號進行2DPSK調(diào)制s_2dpsk=d_sjx.*ht;%生成2DPSK信號%畫出已調(diào)信號2DPSK及其功率譜figure(5

8、)subplot(211)plot(t,s_2dpsk);%畫出2DPSK波形axis(0 10 -1.2 1.2);xlabel('t');ylabel('2DPSK');subplot(212)s_2dpskf,s_2dpsk,df1,f=T2F4(s_2dpsk,ts,df,fs);%求2DPSK信號功率譜plot(f,10*log10(abs(fftshift(s_2dpskf).2/length(f);% 畫出2DPSK信號功率譜axis(-fc-3*B fc+3*B -100 10);xlabel('f');ylabel('2

9、DPSK信號PDF'); %-將2DPSK信號送入信道進行傳輸，先生成信道加性高斯白噪聲噪聲snr_lin=10(snr/10); %換算成倍數(shù)signal_energy=0.5*A2*Ts;%求出接收信號平均能量noise_power=( signal_energy *fs)/(snr_lin*4);%求出噪聲方差（噪聲均值為0）noise_std=sqrt(noise_power);%求出噪聲均方差noise=noise_std.*randn(1,Lt);%以噪聲均方差作為幅度產(chǎn)生高斯白噪聲%-將已調(diào)信號送入信道%畫出信道中的高斯白噪聲及其功率譜figure(6)subplot(2

10、11)plot(t,noise(1:Lt);%畫出噪聲xlabel('t');ylabel('信道噪聲');axis(0 10 -3 3);noisef,noise,df1,f=T2F4(noise,ts,df,fs);%求信道噪聲功率譜subplot(212)plot(f,10*log10(abs(fftshift(noisef).2/length(f);% 畫出信道噪聲功率譜axis(-fs/2-2 fs/2+2 -50 0);xlabel('f');ylabel('信道噪聲PDF');r=s_2dpsk(1:Lt)+nois

11、e(1:Lt);%疊加了噪聲的已調(diào)信號，相當于將已調(diào)信號送入理想信道%畫出加噪后的已調(diào)信號2DPSK及其功率譜figure(7)subplot(211)plot(t,r);%畫出加噪后的已調(diào)信號2DPSKxlabel('t');ylabel('加噪2DPSK信號');axis(0 10 -3 3);rf,r,df1,f=T2F4(r,ts,df,fs);%求加噪后的已調(diào)信號2DPSK功率譜subplot(212)%畫出加噪后已調(diào)信號的功率譜plot(f,10*log10(abs(fftshift(rf).2/length(f);% 畫出已調(diào)信號2DPSK功率譜a

12、xis(-fc-3*B fc+3*B -50 0);xlabel('f');ylabel('加噪2DPSK信號PDF'); %-在接收端準備進行解調(diào)，先通過帶通濾波器%畫出帶通濾波器H,f=bp_f(length(rf),f_start,f_cutoff,df1,fs,1);%經(jīng)過理想帶通濾波器figure(8)subplot(322)plot(f,abs(fftshift(H);% 畫出理想帶通濾波器axis(-fc-3*B fc+3*B -0.2 1.2);xlabel('f');ylabel('理想BPF');DEM = H

13、.*rf; %濾波器輸出的頻譜dem=F2T(DEM,fs);%濾波器的輸出波形dem1=dem(1:Lt);%經(jīng)過理想帶通濾波器后的信號波形及功率譜subplot(323)%經(jīng)過理想帶通濾波器后的信號波形plot(t,dem1);%畫出經(jīng)過理想帶通濾波器后的信號波形axis(0 10 -4 4);xlabel('t');ylabel('理想BPF輸出信號');demf1,dem1,df1,f=T2F4(dem1,ts,df,fs);%求經(jīng)過理想帶通濾波器后信號功率譜subplot(324)plot(f,10*log10(abs(fftshift(demf1).

14、2/length(f);% 畫出經(jīng)過理想帶通濾波器后信號功率譜axis(-fc-3*B fc+3*B -50 0);xlabel('f');ylabel('理想BPF輸出信號PDF');%-進行相干解調(diào)，先和本地載波相乘,即混頻subplot(325);%畫出同頻同相的本地載波plot(t,ht);axis(0 10 -1.2 1.2);xlabel('t');ylabel('本地載波');subplot(326)% 本地載波頻譜htf,ht,df1,f=T2F4(ht,ts,df,fs);plot(f,fftshift(abs(

15、htf);% 畫出載波頻譜axis(-fc-3*B fc+3*B 0 15);xlabel('f');ylabel('本地載波頻譜'); %畫出混頻后的信號及功率譜figure(9) der=dem1(1:Lt).*ht(1:Lt);%和本地載波相乘，即混頻subplot(221)%畫出混頻后的波形plot(t,der);axis(0 10 -2 2);xlabel('t');ylabel('混頻后的信號');derf,der,df1,f=T2F4(der,ts,df,fs);%求混頻后信號的功率譜subplot(222)plot

16、(f,10*log10(abs(fftshift(derf).2/length(f);%畫出混頻后的功率譜axis(-fc-3*B fc+3*B -50 0);xlabel('f');ylabel('混頻后信號的PDF');%-再經(jīng)過低通濾波器%畫出理想低通濾波器LPF,f=lp_f(length(derf),B,df1,fs,1);%求低通濾波器subplot(224) % 畫出理想低通濾波器plot(f,fftshift(abs(LPF);axis(-fc-3*B fc+3*B -0.2 1.2);xlabel('f');ylabel(

17、9;理想LPF');%混頻信號經(jīng)理想低通濾波器后的波形及功率譜DM = LPF.*derf; %理想低通濾波器輸出信號頻譜dm=F2T(DM,fs); %理想低通濾波器的輸出波形 figure(10)subplot(321)plot(t,dm(1:Lt);%畫出經(jīng)過低通濾波器后的解調(diào)出的波形axis(0 10 -1.2 1.2);xlabel('t');ylabel('LPF輸出信號');subplot(322)dmf,dm,df1,f=T2F4(dm,ts,df,fs);%求LPF輸出信號的功率譜plot(f,10*log10(abs(fftshift

18、(dmf).2/length(f);%畫出LPF輸出信號的功率譜axis(-fc-3*B fc+3*B -50 0);xlabel('f');ylabel('LPF輸出信號功率譜');%-最后對LPF輸出信號抽樣判決dm= dm(1:Lt);panjue=zeros(1,N);%建立存儲判決值的矩陣%抽樣判決，規(guī)則：大于等于0判1，小于0判-1for i=1:N; if dm(fc*N_sample*(i-1)+fc*N_sample/2+1)>=0;%抽樣判決時刻 panjue(i)=1; else panjue(i)=-1; endend%-畫出判決出

19、的基帶信號波形，并和調(diào)制信號比較rr=sigexpand(panjue,fc*N_sample);rrt=ones(1,fc*N_sample); % NRZ波形huifu_NRZ=conv(rr,rrt); %觀察2DPSK調(diào)制信號和初步解調(diào)信號波形subplot(323)plot(t,d_sjx(1:Lt);%調(diào)制信號波形axis(0 N -1.2 1.2);xlabel('t');ylabel('調(diào)制信號');subplot(324)plot(t,huifu_NRZ(1:Lt);%解調(diào)信號波形axis(0 N -1.2 1.2);xlabel('t

20、');ylabel('解調(diào)信號');subplot(325)plot(t,d_NRZ(1:Lt);%原始信號波形axis(0 N 0 1.2);xlabel('t');ylabel('原始信號');subplot(326); huifu=mafanbianhuan(huifu_NRZ(1:Lt)+ones(1,Lt)/2,fs);plot(t,huifu);%解調(diào)信號波形axis(0 N 0 1.2);xlabel('t');ylabel('2DPSK解調(diào)信號'); huifu2=;%放置恢復(fù)出來的初始序列

21、值%恢復(fù)原始序列長度for i=0:length(huifu)/16-1 huifu2(i+1)=huifu(i*16+1);end huifu=huifu2; %D/A figure(11);subplot(211);out2=pcm_decode(huifu,2048);t1=0:t0/(length(out2)-1):t0;stem(t1,out2);title('譯碼結(jié)果');axis(0 t0 -3 2.5); y=;%對out2進行擴展，方便畫圖for i=1:length(out2) y=y ones(1,t0/dt1)*out2(i);endt1=0:t0/(l

22、ength(y)-1):t0; dt=0.01;t2=-0.1:dt:0.1;gt=sinc(fs*t2);stt=sigexpand(out2,sdt/dt);out3=conv(stt,gt);subplot(212); t3=0:t0/(length(out3)-1):t0;plot(t3,out3);axis(0 t0 -2.5 2.5);title('通過LPF后的恢復(fù)信號'); %-統(tǒng)計誤碼數(shù)numoferr=sum(abs(panjue-d)/2)/N%計算出錯誤碼元數(shù)functionH,f=bp_f(n,f_start,f_cutoff,df1,fs,p)%帶通

23、濾波器函數(shù) 輸入設(shè)計的濾波器參數(shù)，產(chǎn)生帶通濾波器頻率特性函數(shù)H和頻率向量f%-輸入?yún)?shù)%n 帶通濾波器的輸入信號長度%f_start 通帶起始頻率%f_cutoff 帶通濾波器的截止頻率%df1 頻率分辨率%fs 抽樣頻率%p 濾波器幅度%-輸出(返回)參數(shù)%H 帶通濾波器頻率響應(yīng)%f 頻率向量 %設(shè)計濾波器 n_cutoff = floor(f_cutoff/df1); n_start = floor(f_start/df1); f = 0:df1:df1*(n-1) -fs/2; %頻率向量H = zeros(size(f);H(n_start+1:n_cutoff) = p*ones(1

24、,n_cutoff-n_start);H(length(f) - n_cutoff+1:length(f)-n_start) = p*ones(1,n_cutoff-n_start); function H,f=lp_f(n,f_cutoff,df1,fs,p)%低通濾波器函數(shù) 輸入設(shè)計的濾波器參數(shù)，產(chǎn)生低通濾波器頻率特性函數(shù)H和頻率向量f%-輸入?yún)?shù)%n 低通濾波器的輸入信號長度%f_cutoff 低通濾波器的截止頻率%df1 頻率分辨率%fs 抽樣頻率%p 濾波器幅度%-輸出(返回)參數(shù)%H 低通濾波器頻率響應(yīng)%f 頻率向量n_cutoff = floor(f_cutoff/df1); %

25、設(shè)計濾波器f = 0:df1:df1*(n-1) -fs/2; %頻率向量H = zeros(size(f);H(1:n_cutoff) = p*ones(1,n_cutoff);H(length(f) - n_cutoff+1:length(f) = p*ones(1,n_cutoff);function numoferr,panjue,desingal,t=bpskberr(A,fc,snr,N_sample,N,Ts,d,df)%求誤碼率%-系統(tǒng)仿真參數(shù)% A; %載波振幅% fc; %載波頻率（Hz）% snr; %解調(diào)器輸入信噪比dB% N_sample;%每個碼元的的采樣點數(shù)% N

26、; % 碼元數(shù)% Ts; % 碼元寬度%d;輸入二進制代碼%df:頻率分辨率% -輸出(返回)參數(shù)%numoferr;%誤碼率%panjue 恢復(fù)的二進制代碼1用1表示，0用-1表示%desingal;%恢復(fù)的數(shù)字基帶信號%t;時域采樣時序%-生成調(diào)制信號B=1/Ts;f_start=fc-B;f_cutoff=fc+B;fs=fc*N_sample;%系統(tǒng)采樣頻率ts=Ts/fs; % 系統(tǒng)采樣間隔t=0:ts:N*Ts-ts;Lt=length(t);% 產(chǎn)生二進制信源dd=sigexpand(d+1)/2,fc*N_sample);gt=ones(1,fc*N_sample); % NR

27、Z波形d_NRZ=conv(dd,gt);d_sjx=2*d_NRZ-1;%生成雙極性NRZ信號%-對數(shù)字基帶信號進行2PSK調(diào)制ht=A*sin(2*pi*fc*t);%載波s_2psk=d_sjx(1:Lt).*ht;%生成2PSK信號%-生成高斯白噪聲噪聲snr_lin=10(snr/10); %換算成倍數(shù)signal_energy=0.5*A2*Ts;%求出接收信號平均能量noise_power=( signal_energy *fs)/(snr_lin*4);%求出噪聲方差（噪聲均值為0）noise_std=sqrt(noise_power);%求出噪聲均方差noise=noise_

28、std.*randn(1,Lt);%以噪聲均方差作為幅度產(chǎn)生高斯白噪聲%-將已調(diào)信號送入信道r=s_2psk(1:Lt)+noise(1:Lt);%疊加了噪聲的已調(diào)信號，相當于將已調(diào)信號送入理想信道rf,r,df1,f=T2F4(r,ts,df,fs);%-在接收端先通過帶通濾波器H,f=bp_f(length(rf),f_start,f_cutoff,df1,fs,1);DEM = H.*rf; %濾波器輸出的功率譜dem=F2T(DEM,fs);%濾波器的輸出信號波形 %進行相干解調(diào)，先和本地載波相乘,即混頻%和本地載波相乘der=dem(1:Lt).*ht;%混頻derf,der,df1

29、,f=T2F4(der,ts,df,fs);%求混頻后信號的功率譜%-再經(jīng)過低通濾波器LPF,f=lp_f(length(derf),B,df1,fs,1);DM = LPF.*derf; %理想低通濾波器輸出的功率譜dm=F2T(DM,fs);%理想低通濾波器的輸出信號%-最后對LPF輸出信號抽樣判決panjue=zeros(1,N);%建立存儲判決值的矩陣%抽樣判決，規(guī)則：大于等于0判1，小于0判-1for i=1:N if dm(fc*N_sample*(i-1)+fc*N_sample/2+1)>=0;%抽樣判決時刻 panjue(i)=1; else panjue(i)=-1;

30、 endend%-生成判決出的基帶信號波形dd1=sigexpand(panjue,fc*N_sample);gt1=ones(1,fc*N_sample); % NRZ波形desinga=conv(dd1,gt1);desingal=desinga(1:Lt);%-統(tǒng)計誤碼數(shù)numoferr=sum(abs(panjue-d)/2)/N;%計算出錯誤碼元數(shù)function y=dpskberr(d)%-不同信噪比下的誤碼率 A=1; %載波振幅fc=2; %載波頻率（Hz）SNRindB1=-5:1:12; % 信噪比取值向量，dB為單位SNRindB2=-5:0.2:12; % 信噪比取值

31、向量，dB為單位%N_sample=100;%每個碼元的的采樣點數(shù)N_sample=50;N=2000; % 碼元數(shù)d=sign(rand(1,N)-0.5+eps);%產(chǎn)生雙極性二進制碼元 Ts=1; % 碼元寬度df=0.01;simu_err_prb=zeros(1,length(SNRindB1);%理論誤碼率for i=1:length(SNRindB2) % 計算信噪比值 SNR=10(SNRindB2(i)/10); % 計算普通接收機的理論誤碼率。 theo_err_prb(i)= 1-erf(sqrt(SNR)2; %一般蒙特卡洛仿真是針對最佳接收機的，對于來自同一信道的接收

32、信號， %信道高斯噪聲的功率譜密度是不變的，最佳接收機性能之所以比一般接收機好， %是由于其輸入信噪比是一般接收機的2倍（一般接收機帶通濾波器帶寬按2/Ts計算）； %故信道高斯白噪聲的單邊帶功率譜密度即可按最佳接收機設(shè)計，也可按普通接收機設(shè)計。 %互補誤差函數(shù) theo_err_prb1(i)= 1-erf(sqrt(2*SNR)2;%最佳接收機理論誤碼率曲線end%計算普通接收機實際誤碼率for i=1:length(SNRindB1) % d 矩陣為模擬信號經(jīng)過抽樣編碼，碼變換后的矩陣 numoferr,panjue,desingal,t=bpskberr(A,fc,SNRindB1(i

33、),N_sample,N,Ts,d,df); simu_err_prb(i)=numoferr;end% 誤碼率曲線圖：估計值和理論值曲線對比圖figure(12);semilogy(SNRindB2,theo_err_prb,SNRindB1,simu_err_prb,'o',SNRindB2,theo_err_prb1,'*');axis(-5 12 0.00000001 1);xlabel('SNR in dB');ylabel('Prb of Err');legend('普通接收機理論誤碼率','普

34、通接收機實際誤碼率','最佳接收機理論誤碼率');function m=F2T(M,fs)%-輸入?yún)?shù)%M：信號的頻譜%fs:系統(tǒng)采樣頻率%-輸出(返回)參數(shù)%m:傅里葉逆變換后的信號，注意其長度為2的整數(shù)次冪，利用其畫波形時，要注意選取m的一部分，選取長度和所給時間序列t的長度要一致，plot(t,m(1:length(t)，否則會出錯。m = real(ifft(M)*fs; function M,m,df=fftseq(m,ts,df)%各參數(shù)含義與子函數(shù)T2F中的完全相同，完成fs = 1/ts;if nargin =2 n1 =0;else n1 = fs/d

35、f;end n2 = length(m);n = 2(max(nextpow2(n1),nextpow2(n2);M = fft(m,n);m = m,zeros(1,n-n2);df = fs/n;%碼變換function d_NRZ2=mabianhuan(d_NRZ1,fs)b=1;max=length(d_NRZ1)/fs;for m=1:max n=1+(m-1)*fs; d_NRZ2(n:n+fs-1)=mod2(d_NRZ1(n),b)*ones(1,fs); b=d_NRZ2(n);end%碼反變換function yuanshi=mafanbianhuan(jietiao,f

36、s)b=1;max=length(jietiao)/fs;for m=1:max n=1+(m-1)*fs; yuanshi(n:n+fs-1)=mod2(jietiao(n),b)*ones(1,fs); b=jietiao(n);end%模2加運算function y=mod2(a,b)if a+b=2 y=0;else y=a+b;end%折線譯碼函數(shù)：%輸入in為8位二進制碼，（v，+v）為量化區(qū)間function out=pcm_decode(in,v)n=length(in);in=reshape(in',8,n/8)'slot(1)=0;slot(2)=16;sl

37、ot(3)=32;slot(4)=64;slot(5)=128;slot(6)=256;slot(7)=512;slot(8)=1024;step(1)=1;step(2)=1;step(3)=2;step(4)=4;step(5)=8;step(6)=16;step(7)=32;step(8)=64;for i=1:n/8 ss=2*in(i,1)-1; tmp=in(i,2)*4+in(i,3)*2+in(i,4)+1; st=slot(tmp);dt=(in(i,5)*8+in(i,6)*4+in(i,7)*2+in(i,8)*step(tmp)+0.5*step(tmp); out(i

38、)=ss*(st+dt)/2048*v;end function out=pcm_encode(x)n=length(x);for i=1:n%編寫極性碼 if x(i)>0 out(i,1)=1; else out(i,1)=0; end%編寫段落碼并計算量化間隔和量化起始電平 if 0<=abs(x(i)&abs(x(i)<16%段落碼out(i,2)=0;out(i,3)=0;out(i,4)=0;%量化間隔step=1;%起始電平st=0; elseif 16<=abs(x(i)&abs(x(i)<32out(i,2)=0;out(i,3)

39、=0;out(i,4)=1;step=1;st=16; elseif 32<=abs(x(i)&abs(x(i)<64out(i,2)=0;out(i,3)=1;out(i,4)=0;step=2;st=32; elseif 64<=abs(x(i)&abs(x(i)<128out(i,2)=0;out(i,3)=1;out(i,4)=1;step=4;st=64; elseif 128<=abs(x(i)&abs(x(i)<256out(i,2)=1;out(i,3)=0;out(i,4)=0;step=8;st=128; elseif 256<=abs(x(i)&abs(x(i)<512out(i,2)=1;out(i,3)=0;out(i,4)=1;step=16;st=256; elseif 512<=abs(x(i)&abs(x(i)<1024out(i,2)=1;out(i,3)=1;out(i,4)=0;step=32;st=512;el