2DPSK调制与解调matlab.doc

2DPSK 调制与解调
%---------------------------------------------------
%>>>>>>>>>>>>>>>>>>参数初始化>>>>>>>>>>>>>>>>>>>>>
%---------------------------------------------------
fs = 3600000;%采样频率为36000赫兹
Time_Hold_On = 1/1200;%一个时钟周期为1200分之1，对应比特率为1200bps
Num_Unit = fs * Time_Hold_On;%一个时钟周期内的采样点个数
High_Level = ones ( 1, Num_Unit );%高电平（全1序列）
Low_Level = zeros ( 1, Num_Unit );%低电平（全0序列）
w = 1800;%载波角频率1800Hz
A = 1;%载波幅值
%---------------------------------------------------
%>>>>>>>>>>>>>>>>>>信号初始化>>>>>>>>>>>>>>>
%---------------------------------------------------
Sign_Set = [0,1,1,0,1,0,0,1];%原始序列
Lenth_Of_Sign = length ( Sign_Set );%原始序列长度
Sign_Sett = ones(1,Lenth_Of_Sign+1);%差分变换后的序列，初始化为长度为原始序列长度+1的全1序列（第一个码元为1）
sign_orign = zeros ( 1, Num_Unit * (Lenth_Of_Sign+1) );%初始化基带信号为全0序列
sign_result = zeros ( 1, Num_Unit * (Lenth_Of_Sign+1) );%初始化接收到的基带信号为全0序列
st = zeros ( 1, Num_Unit *( Lenth_Of_Sign+1) );%初始化调制后的信号为全0序列
t = 0 : 1/fs : Time_Hold_On * (Lenth_Of_Sign +1)- 1/fs;%信号采样时间点
result=zeros(1,Lenth_Of_Sign+1);%初始化接收到的序列
resultt=zeros(1,Lenth_Of_Sign);%初始化差分解调后的序列

%---------------------------------------------------
%>>>>>>>>>>>求差分编码>>>>>>>>>>>>
%-----------------

for I = 2 : Lenth_Of_Sign+1 %差分变换后的序列第一个值为1，从第2个开始计算
Sign_Sett(I)= xor(Sign_Sett(I-1),Sign_Set(I-1));%用异或运算求差分码
end

%---------------------------------------------------
%>>>>>>>>>>>产生基带信号>>>>>>>>>>>>
%---------------------------------------------------

for I = 1 : Lenth_Of_Sign+1 %考虑差分变换后序列中每一个值
if Sign_Sett(I) == 1
sign_orign( (I-1)*Num_Unit + 1 : I*Num_Unit) = High_Level; %序列值为1，基带信号为高电平
else
sign_orign( (I-1)*Num_Unit + 1 : I*Num_Unit) = Low_Level; %序列值为0，基带信号为低电平
end
end

%------------------------------------------------
%>>>>>>>>>>>>>>>>>>调制部分>>>>>>>>>>>>>>>>>>
%---------------------------------------------------
for I = 1 : Lenth_Of_Sign+1 %考虑差分变换后序列中每一个值
if Sign_Sett(I) == 1
st( (I-1)*Num_Unit + 1 : I*Num_Unit) = A * cos ( 2 * pi * w * t( (I-1)*Num_Unit + 1 : I*Num_Unit ) + ( pi / 2 ) );%序列值为1，相位调制为π/2
else
st( (I-1)*Num_Unit + 1 : I*Num_Unit) = A * cos ( 2 * pi * w * t( (I-1)*Num_Unit + 1 : I*Num_Unit ) );%序列值为0，相位调制为0
end
end

figure
subplot ( 2,1,1 )
plot(t, sign_orign);
axis( [ 0 , Time_Hold_On *( Lenth_Of_Sign + 2), - (A / 2), A + (A / 2) ] );
title ( '原始信号' );
grid %画出基带信号

subplot ( 2, 1, 2 );
plot ( t, st );
axis( [ 0 , Time_Hold_On *( Lenth_Of_Sign + 2), - 3*(A / 2), 3*(A / 2) ] );
title ( '调制后的信号' );
grid %画出相位调制后的信号

%---------------------------------------------------
%>>>>>>>>>>>>>>>>>>相干解调>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%--------------------------
dt = st .* cos ( 2 * pi * w * t ); %相干相乘

figure
plot ( t, dt );
axis( [ 0 , Time_Hold_On *( Lenth_Of_Sign + 2), - 3*(A / 2), 3*(A / 2) ] );
title ( '相干相乘后的波形' );
grid

%---------------------------------------------------
%>>>>>>>>>>>>>>>>>>>低通滤波部分>>>>>>>>>>>>>>>>>>>>
%---------------------------------------------------
[N,Wn] = buttord( 2*pi*1500, 2*pi*3400,3,25,'s'); %临界频率采用角频率表示，计算低通滤波器参数
[b,a]=butter(N,Wn,'s'); %产生N阶低通巴特沃斯滤波器
[bz,az]=impinvar(b,a,fs); %映射为数字的

dt = filter(bz,az,dt); %将相干相乘后的信号进行滤波
figure
plot ( t, dt );
axis( [ 0 , Time_Hold_On *( Lenth_Of_Sign + 2), - 3*(A / 2), 3*(A / 2) ] );
title ( '低通滤波后的波形' );
grid

%---------------------------------------------------
%>>>>>>>>>>>>>抽样判决 & 逆码变换部分>>>>>>>>>>>>>>>
%---------------------------------------------------
for I = 1 : Lenth_Of_Sign+1
if dt((2*I-1)*Num_Unit/2) < 0.25 在时钟周期中间采样，由于相干相乘后信号幅值变为1/2，所以判决门限为0.25
sign_result( (I-1)*Num_Unit + 1 : I*Num_Unit) = High_Level;
else
sign_result( (I-1)*Num_Unit + 1 : I*Num_Unit) = Low_Level;
end
%由于相干相乘后得到的信号变为原来的负数，所以进行逆码变换
end
figure
plot ( t, sign_result );
axis( [ 0 , Time_Hold_On *( Lenth_Of_Sign + 2), - 3A / 2), 3*(A / 2) ] );
title ( '逆码变换后的波形' );
grid

%---------------------------------------------------
%>>>>>>>>>>>序列生成与差分还原>>>>>>>>>>>>
%---------------------------------------------------
for I=1: Lenth_Of_Sign+1
result(I)=sign_result(I* Num_Unit)%将接收到的信号生成（差分）序列
end

for I = 1 : Lenth_Of_Sign
resultt(I)= xor(result(I),result(I+1)) %将接收到的差分序列还原为原序列
end

window=boxcar(length(st)); %矩形窗  
nfft=1024;  
[Pxx,f]=periodogram(st,window,nfft,fs); %求功率谱密度
plot(f,10*log10(Pxx));