gmsk.m

无线通信数字信号GMSK调制解调并画出星座图

%-------------------------------
%    Matlab仿真GMSK调制解调
%          
%------------------------------

N=1000;                             %码元数
S=16;                               %每个码元采样数
Rb=1;                               %码元数率1kbps
Tb=1/Rb;
t=Tb/S;                             %采样时间间隔
fc=4;                               %载波4kHz
BT = 1;                   	    %高斯滤波器BT值，通过改变此值改变功率谱
n = 40;                    			%高斯滤波器参数
o = 16;                    			%高斯滤波器参数，采样频率

num =Rb* N+1;                       %数据长度   

gmsk = zeros(1,N*S+1); 		        %gmsk总长度

%产生双极性基带信号
ck=randint(num,1);      
for k0=1:num
    if ck(k0)==0
        ck(k0)=-1;
    end
end
ck(1) = -1;

Qkt = zeros(1,num);              %正交数据初始化
Ikt = Qkt;                       %同相数据初始化

%串并转换
Ikt(1:2:num-2)=ck(2:2:num-1);    
lkt(2:2:num-1)=ck(2:2:num-1);

Qkt(2:2:num-1)=ck(3:2:num);
Qkt(3:2:num)=ck(3:2:num);

Ikt(num) = 0;

%升采样
Ik = zeros(1,num*S);    		
Qk = Ik;
for k1=1:num
    for k2=1:S
        Ik(S*(k1-1)+k2)=Ikt(k1);
        Qk(S*(k1-1)+k2)=Qkt(k1);
    end
end

 %通过高斯滤波器
gIk = conv(Ik,gaussfir(BT,n,o));
gQk = conv(Qk,gaussfir(BT,n,o));
lgh = length(gIk);

gI=gIk.*cos(pi*t*(1:lgh)*Rb/2).*cos(fc*2*pi*t*(1:lgh));  %gmsk同相分量
gQ=gQk.*sin(pi*t*(1:lgh)*Rb/2).*sin(fc*2*pi*t*(1:lgh));  %gmsk正交分量
 
gmsk= gI + gQ ;  %gmsk时域信号波形

y = fft(gmsk);                  	%生成gmsk信号功率谱
Sgmsk = 10 * log10(y.*conj(y));
a=length(Sgmsk);

%显示gmsk信号功率谱
figure(1) 
plot((0:a/2-fc*a/S)/a*S,Sgmsk(fc*a/S:a/2));
title('gmsk信号功率谱波形');
xlabel('(f-fc)T'), ylabel('dB');


%瑞利多径衰落信道
r=gmsk;
for k=1:50          %假设有50个信道
    b=raylrnd(1);   %瑞利分布的乘性噪声
    delay=floor(S*rand(1)/10);  %延时
    rm=b*[zeros(1,delay),gmsk(1:lgh-delay)];
    r=r+rm;
end

%双边带宽为2Rb的LPF
w=Rb*t*lgh;
LPF=[ones(1,w),zeros(1,lgh-2*w),ones(1,w)];

%同向分量提取
ri=r.*cos(2*pi*fc*t*(1:lgh));
ri=real(ifft(fft(ri,lgh).*LPF));

%正交分量提取
rq=r.*sin(2*pi*fc*t*(1:lgh));
rq=real(ifft(fft(rq,lgh).*LPF));

%取样判决
rik=ri.*sign(cos(pi/2/Tb*t*(1:lgh)));
rqk=rq.*sign(sin(pi/2/Tb*t*(1:lgh)));

result=zeros(1,lgh);

for k=1:N/2                 %并串转换
   result((2*k-2)*S+1:(2*k-1)*S)=sign(rqk((2*k-1)*S));
   result((2*k-1)*S+1:2*k*S)    =sign(rik(2*k*S));
end



figure(2)
plot(t:t:10*S*t,result(1:10*S));
axis([0,10*Tb,-1.5,1.5]);
title('接收到的时域波形');
xlabel('t(ms)');

ber = 0;
for k3 = 1:N-1				%计算误码率
    if result(k3) ~= ck(k3)
        ber = ber + 1;
    end
end
ber