实验一 基于matlab的GMSK综合实训  下载本文

F_sample=64; %每载波采样64个点 B_num=8; %基带信号为8个码元

B_sample=F_sample*Fc*Tb; %每基带码元采样点数_sample=Tb/Dt

Dt=1/Fc/F_sample; %采样间隔 t=0:Dt:B_num*Tb-Dt; %仿真时间 T=Dt*length(t); %仿真时间值 Ak=[0 0 1 0 1 0 1 0]; %产生8个基带信号 Ak=2*Ak-1;

gt=ones(1,B_sample); %每码元对应的载波信号 Akk=sigexpand(Ak,B_sample); %码元扩展 temp=conv(Akk,gt); %码元扩展 Akk=temp(1:length(Akk)); %码元扩展

tt=-2.5*Tb:Dt:2.5*Tb-Dt;

%g(t)=Q[2*pi*Bb*(t-Tb/2)/sqrt(log(2))]-Q[2*pi*Bb*(t+Tb/2)/sqrt(log(2))]; %Q(t)=erfc(t/sqrt(2))/2;

gausst=erfc(2*pi*Bb*(tt-Tb/2)/sqrt(log(2))/sqrt(2))/2-erfc(2*pi*Bb*(tt+Tb/2)/sqrt(log(2))/sqrt(2))/2;

J_g=zeros(1,length(gausst)); %使J_g 的长度和Gausst的一样 for i=1:length(gausst) if i==1

J_g(i)=gausst(i)*Dt; else

J_g(i)=J_g(i-1)+gausst(i)*Dt; end; end;

J_g=J_g/2/Tb; %计算相位Alpha

Alpha=zeros(1,length(Akk)); k=1; L=0;

for j=1:B_sample

J_Alpha=Ak(k+2)*J_g(j);

11

Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2; end; k=2; L=0;

for j=1:B_sample

J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample); Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2; end; k=3; L=0;

for j=1:B_sample

J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample);

Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2; end; k=4; L=0;

for j=1:B_sample

J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j+3*B_sample);

Alpha((k-1)*B_sample+j)=pi*J_Alpha+L*pi/2; end; L=0;

for k=5:B_num-2 if k==5 L=0; else

L=L+Ak(k-3); end;

12

for j=1:B_sample

J_Alpha=Ak(k+2)*J_g(j)+Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j+3*B_sample)+Ak(k-2)*J_g(j+4*B_sample); Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2; end; end;

%B_num-1; k=B_num-1; L=L+Ak(k-3); for j=1:B_sample

J_Alpha=Ak(k+1)*J_g(j+B_sample)+Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j +3*B_sample)+Ak(k-2)*J_g(j+4*B_sample);

Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2; end; %B_num; k=B_num; L=L+Ak(k-3); for j=1:B_sample

J_Alpha=Ak(k)*J_g(j+2*B_sample)+Ak(k-1)*J_g(j+3*B_sample)+Ak(k-2)*J_g(j+4*B_sample);

Alpha((k-1)*B_sample+j)=pi*J_Alpha+mod(L,4)*pi/2; end;

S_Gmsk=cos(2*pi*Fc*t+Alpha); subplot(311) plot(t/Tb,Akk); axis([0 8 -1.5 1.5]); title('基带波形');

subplot(312) plot(t/Tb,Alpha*2/pi);

13

axis([0 8 min(Alpha*2/pi)-1 max(Alpha*2/pi)+1]); title('相位波形');

subplot(313) plot(t/Tb,S_Gmsk); axis([0 8 -1.5 1.5]); title('GMSK波形'); %解调 for n=1:512; if n<=B_sample Alpha1(n)=0;

else Alpha1(n)=Alpha(n-B_sample); end end

a=[0 1 1 1 1 1 1 1 ]

ak=sigexpand(a,B_sample); %码元扩展 temp=conv(ak,gt); %码元扩展 ak=temp(1:length(ak));

S_Gmsk1=cos(2*pi*Fc*(t-Tb)+Alpha1+pi/2).*ak; %延迟1bt,移相pi/2 figure subplot(311) plot(t/Tb,S_Gmsk1); axis([0 8 -1.5 1.5]);

title('延迟1bt,移相pi/2GMSK波形');

xt=S_Gmsk1.*S_Gmsk; x=0; subplot(312) plot(t/Tb,xt,t/Tb,x,'r:'); axis([0 8 -1.5 1.5]); title('相乘后波形'); %低通滤波 Fs=10000; rp=3;rs=50;

14

wp=2*pi*50;ws=2*pi*800; [n,wn]=buttord(wp,ws,rp,rs,'s') [z,p,k]=buttap(n); [bp,ap]=zp2tf(z,p,k); [bs,as]=lp2lp(bp,ap,wn); [b,a]=bilinear(bs,as,Fs) y=filter(b,a,xt); subplot(313) plot(t/Tb,y,t/Tb,x,'r:'); axis([0 8 -1.5 1.5]);

title('经过低通滤波器后波形'); for i=1:8

if y(i*B_sample)>0 bt(i)=1 else bt(i)=0 end end bt=2*bt-1;

btt=sigexpand(bt,B_sample); temp1=conv(btt,gt); btt=temp1(1:length(btt)); figure subplot(311) plot(bt) title('抽样值'); axis([0 8 -1.5 1.5]); subplot(312) plot(t/Tb,Akk); axis([0 8 -1.5 1.5]); title('原基带波形'); subplot(313) plot(t/Tb,btt); axis([0 8 -1.5 1.5]);

%码元扩展 %码元扩展%码元扩展

15

title('解调后波形');

[f F_Gmsk]=T2F(t,S_Gmsk);

plot(f/Fc,10*log10(abs(F_Gmsk).^2/length(F_Gmsk)),'k'); %axis([-10 10 -80 -10]); xlabel('f');

title('GMSK功率谱密度');

以下代码是上面的sigexpand函数的m文件的代码 function [out]=sigexpand(d,M) N=length(d); out=zeros(M,N); out(1,:)=d;

out=reshape(out,1,M*N);

下面的是求功率谱密度的T2F函数的m文件代码 function [f,sf]=T2F(t,st)

%这个子程序主要是运用FFT对一个信号做傅立叶变换 %输入是时间和信号的变量,长度必须大于2 %输出是频谱 dt=t(2)-t(1); T=t(end); df=1/T; N=length(st);

f=-N/2*df:df:N/2*df-df; sf=fft(st); sf=T/N*fftshift(sf);

16