📄 gmsk_msk.m
字号:
% ===============================================================
% Creat by LuyoMig
% 2007.04.04
% ===============================================================
clc;
clear all;
close all;
Ts=1/16000; % 基带信号周期为1/16000s,即为16KHz
Tb=1/32000; % 输入信号周期为Ts/2=1/32000s,即32KHz
BbTb=0.3; % 取BbTb为0.3
Bb=BbTb/Tb; % 3dB带宽
Fc=32000; % 载波频率为32KHz
F_sample=64; % 每载波采样64个点
Ak=sign(randn(1,100)); % 随机产生100个基带信号
% Ak = [0 1 1 1 1 0 1 1 1 1 1 0 1 0 1 0 1 1 1 1 1 1 1 1 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1];
% Ak = [1 0 1 1 1 0 1 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 1 0 1 1 0 1 1 0 1 1 0 1 1];
% Ak=[0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1];
% Ak=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
% Ak=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
% Ak=[0 0 0 0 0]
B_num=length(Ak) % 基带信号为B_num个码元
B_sample=F_sample*Fc*Tb % 每基带码元采样点数B_sample=Tb/Dt
Dt=1/Fc/F_sample % 采样间隔
t=0:Dt:B_num*Tb-Dt; % 仿真时间
T=Dt*length(t); % 仿真时间值
Ak=1-2*Ak; % 双极性变换
gt=ones(1,B_sample); % 每码元对应的载波信号,不归零码
Akk=sigexpand(Ak,B_sample); % 码元扩展
Akk;
L_Akk=length(Akk)
temp=conv(Akk,gt); % 码元扩展
Akk=temp(1:L_Akk); % 码元扩展
% ===============================================================
% 以下部分计算 MSK 调制出来的波形
% ===============================================================
Tbt=0:Dt:Tb-Dt; % 每Tb时间
ThetaL=0;
Theta=zeros(1,L_Akk);
for k=1:B_num % 计算每个采样点的相位 Theta
if k==1
ThetaL=0;
else
ThetaL=Theta((k-1)*B_sample) % 调制完k个码元时的相位
end
for j=1:B_sample % pi/2*Ak*Tbt/Tb 计算当前相位
Theta((k-1)*B_sample+j)=pi*Ak(k)*Tbt(j)/2/Tb+ThetaL; % 当前相位等于前面所有码元相位的累加
end
end;
Theta;
S_Msk=cos(2*pi*Fc*t+Theta); % 调制到载波上的最后MSK波形
S_Msk_I=cos(Theta); % 基带I路
S_Msk_Q=-sin(Theta); % 基带Q路
[f F_Msk]=T2F(t,S_Msk);
% ===============================================================
% 根据公式计算g(t)
% ===============================================================
tt=-2.5*Tb:Dt:2.5*Tb-Dt; % 把高斯滤波器的响应长度截短长度为5
%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; % 计算g(t)
J_g=zeros(1,length(gausst)); % 使J_g 的长度和Gausst的一样，方便后面进行积分
for i=1:length(gausst) % 积分求和 积分的最后结果等于 1/2
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;
L_J_g=length(J_g)
% ===============================================================
% 以下部分计算 GMSK公式S(t)=cos(2*pi*Fc*t+alpha(t))中的alpha(t)
% ===============================================================
Alpha=zeros(1,length(Akk));
k=1;
L=0; % 定义初始相位都为"0"
for j=1:B_sample
J_Alpha=Ak(k+2)*J_g(j);
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;
% fid=fopen('out.txt','w');
L=0;
for k=5:B_num-2
if k==5
L=0;
else
L=L+Ak(k-3);
end;
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+L*pi/2;
% fprintf(fid,'%f ',Alpha((k-1)*B_sample+j));
% fprintf(fid,'%f ',Alpha((k-1)*B_sample+j)*65535);
% fprintf(fid,'%f ',Alpha((k-1)*B_sample+j)*10^6);
% fprintf(fid,'%d \r\n',L);
end;
end;
% fclose(fid);
%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+L*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+L*pi/2;
end;
% Alpha(7*B_sample)
% Alpha(7*B_sample)*10^6
% ===============================================================
% 以下部分计算 GMSK 最后调制出来的波形
% ===============================================================
S_Gmsk=cos(2*pi*Fc*t+Alpha);
%---------------------------------------------
S_Gmsk_I=cos(Alpha); % 基带I路
S_Gmsk_Q=-sin(Alpha); % 基带Q路
%---------------------------------------------
% ===============================================================
% 以下部分计算 GMSK 最后调制出来波形的频域图
% ===============================================================
[f F_Gmsk]=T2F(t,S_Gmsk);
[f_I F_Gmsk_I]=T2F(t,S_Gmsk_I);
[f_Q F_Gmsk_Q]=T2F(t,S_Gmsk_Q);
figure(1);
subplot(2,1,1);
plot(tt/Tb,gausst);
% plot(tt,gausst);
axis([-2.5 2.5 0 1]);
title('Gaussian');
subplot(2,1,2);
plot(tt/Tb,J_g);
axis([-2.5 2.5 0 0.5]);
title('积分后的值');
%
%
Show_Num=B_num; % 显示码元数
Show_Time=Show_Num*Tb; % 显示码元数
figure(2); % 基带信号
subplot(421)
plot(t/Tb,Akk);
axis([0 Show_Num -1.5 1.5]);
% plot(Akk);
% axis([0 3300 -1.5 1.5]);
title('基带波形');
subplot(423) % MSK相位图
plot(t/Tb,Theta*2/pi);
axis([0 Show_Num min(Theta*2/pi)-1 max(Theta*2/pi)+1]);
title('MSK相位');
grid on;
subplot(425) % MSK调制波形
plot(t/Tb,S_Msk);
axis([0 Show_Num -1.5 1.5]);
title('MSK波形');
subplot(427)
plot(f/Fc,10*log10(abs(F_Msk).^2/length(F_Msk)));
axis([-15 15 -200 -100]);
% plot(f/Fc,10*log10(abs(F_Msk).^2/Tb));
% axis([-50 50 -200 -10]);
xlabel('f');
ylabel('MSK功率谱密度(dB/Hz)');
title('MSK功率谱密度');
subplot(422)
plot(t/Tb,Akk,'k');
axis([0 Show_Num -1.5 1.5]);
title('基带波形');
subplot(424)
plot(t/Tb,Alpha*2/pi,'k');
axis([0 Show_Num min(Alpha*2/pi)-1 max(Alpha*2/pi)+1]);
title('GMSK相位波形');
grid on;
subplot(426)
plot(t/Tb,S_Gmsk,'k');
% plot(t,S_Gmsk);
axis([0 Show_Num -1.5 1.5]);
title('GMSK波形');
subplot(428)
plot(f/Fc,10*log10(abs(F_Gmsk).^2/length(F_Gmsk)),'k');
axis([-15 15 -200 -100]);
% plot(f/Fc,10*log10(abs(F_Gmsk).^2/Tb));
% axis([-50 50 -200 -10]);
xlabel('f');
title('GMSK功率谱密度');
figure(3)
plot(t/Tb,Alpha*2/pi,'k');
axis([0 Show_Num min(Alpha*2/pi)-1 max(Alpha*2/pi)+1]);
title('GMSK相位波形');
grid on;
% figure(5)
% plot(f/Fc,F_Gmsk,'k');
% % axis([-15 15 -200 -100]);
% % plot(f/Fc,10*log10(abs(F_Gmsk).^2/Tb));
% % axis([-50 50 -200 -10]);
% xlabel('f');
% title('GMSK频域图');
%********************************IQ输出2007.06.08***************************
% %---------------数据比较----------------
% load ('I_path_2.mat','-mat');
% load ('Q_path_2.mat','-mat');
% %---------------------------------------
figure(4)
subplot(211)
plot(t/Tb,S_Gmsk_I);
axis([0 Show_Num+1 -1.5 1.5]);
title('GMSK波形 I路');
grid on;
% hold on; % 比较数据 I路
% plot(t/Tb,I_path,'r');
% axis([0 Show_Num+1 -1.5 1.5]);
% legend('理论结果','实际结果');
% grid on;
subplot(212)
plot(t/Tb,S_Gmsk_Q);
axis([0 Show_Num+1 -1.5 1.5]);
title('GMSK波形 Q路');
grid on;
% hold on; % 比较数据 Q路
% plot(t/Tb,Q_path,'r');
% axis([0 Show_Num+1 -1.5 1.5]);
% legend('理论结果','实际结果');
% grid on;
% subplot(223)
% plot(t/Tb,S_Msk_I,'k');
% axis([0 Show_Num -1.5 1.5]);
% title('MSK波形 I路');
%
% subplot(224)
% plot(t/Tb,S_Msk_Q,'k');
% axis([0 Show_Num -1.5 1.5]);
% title('MSK波形 Q路');
% figure
% plot(f_I/Fc,10*log10(abs(F_Gmsk_I).^2/length(F_Gmsk_I)));
% axis([-15 15 -200 -100]);
% xlabel('f');
% ylabel('MSK功率谱密度(dB/Hz)');
% title('MSK功率谱密度');
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -