dbpsk_modem.m

dbpsk信号的调制解调matlab代码实现的一个例子

%解调采用1bit差分延迟检测
%调制时差分编码，这样用1bit延迟差分检测的结果，直接就是原码


clear;
close all;

transmit_data=[0,0,1,0,0,1,0,0,1,0,1,1,0,0,1,1,1,1,0,1,0,1,0,1,0,0,1,0,1,1,0,0,1,0,1,0,1,1,1,0,0,1,0,0,0,0,0,1,1,0,0,1,1,0,0,0,1,0,0,0,1,0,1,0,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,0,1,1,1,1,1,0,0,0,0,1,1,1,0,0,1,0,1,1];
N=length(transmit_data);%数据长度


Tb=1;%码元周期
Fb=1/Tb;

Fc=Fb;%载波频率等于码元速率

Q=32;%一个码元间隔tb里的采样点数
Fs=Q*Fb;
Ts=1/Fs;


%平方根升余弦滤波器冲击响应h_rrc
alpha=0.35;
h_rrc=rcosine(Fb,Fs,'fir/sqrt',alpha,3);%默认的是6个Tb长度
h_rrc=h_rrc(17:177);%截取中间161个点，5个Tb
h_rrc=h_rrc/sum(h_rrc);

%差分编码
dif(1)=0;
for i=1:N;
    if transmit_data(i)==dif(i);
        dif(i+1)=0;
    else dif(i+1)=1;
    end
end
%转换成双极性码
for i=1:N;
    if dif(i)==0;
       dif(i)=-1;
    end
end

%上采样，即插入0点
g=upsample(dif,Q);

%通过根升余弦滤波器
g_low=convn(g,h_rrc,'same');
g_low=g_low/max(g_low);

figure;
plot(g_low);
title('g_low');

%调制
n=0:length(g_low)-1;
dbpsk=g_low.*cos(2*pi*Fc*n*Ts);

figure;
plot(dbpsk);
title('dbpsk');

% %加进噪声
% N1=length(dbpsk);
% dbpsk_spectrum=fft(dbpsk);
% dbpsk_power=(norm(dbpsk)^2)/N1;
% SNRindB=1;%输出信噪比
% SNR=10^(SNRindB/10);
% noise_power=dbpsk_power/(SNR);
% noise=sqrt(noise_power)*randn(1,N1);
% dbpsk_noise=dbpsk+noise;

%截取一段
dbpsk_trans=dbpsk(11:end-7);

%解调
for i=Q+1:length(dbpsk_trans);
    product(i-Q)=dbpsk_trans(i)*dbpsk_trans(i-Q);%原信号和延迟1bit后的信号相乘
end

h_low=fir1(50,0.01,'low');%低通滤波器
product_low=convn(product,h_low,'same');
figure;
plot(product);
title('product');
figure;
plot(product_low);
title('product_low');

%变成双极性不归零码
for i=1:length(product_low);
    if product_low(i)<0;%注意，这里是反的，小于零的才判为1！
        data_reshape(i)=1;
    else data_reshape(i)=-1;
    end
end

swichtime_counter=56;%码元转换时刻计数器初值
delta=8;%delta为Tb/4
ee=5;%vd的容限
i2=1;%为了观察vd的变化，将vd设成数组。在实用中，不必如此，将vd设成一个数值就可以了
n=1;
for i=1:length(data_reshape)-32;
    swichtime_counter=swichtime_counter-1;%码元转换时刻计数器
    if swichtime_counter==0;%码元转换时刻计数器为0时，码元转换时刻到来，开始计算早积分和迟积分
        %计算早积分
        acc=0;
        for j=1:32;
            acc=acc+data_reshape(i-(Q-delta)+j);
        end
        early_integral(i2)=abs(acc);%取绝对值
        
        %计算迟积分
        acc=0;
        for j=1:32;
            acc=acc+data_reshape(i-delta+j);
        end
        late_integral(i2)=abs(acc);%取绝对值
        
        vd(i2)=early_integral(i2)-late_integral(i2);%得到误差信号vd
        %根据误差信号vd来改变码元转换时刻计数器值
        if vd(i2)>ee;
            swichtime_counter=33;
        elseif vd(i2)<-ee;
            swichtime_counter=31;
        else swichtime_counter=32;
        end
        swichtime(i)=2;%为了观察取样时刻，设立的一个临时量。
        i2=i2+1;
        
        %采样
        receive_data(n)=data_reshape(i-16);
        if receive_data(n)==-1;
            receive_data(n)=0;
        end
        n=n+1;
    end
end


figure;
stem(data_reshape);
title('data_reshape');
hold on;
stem(swichtime,'r','.');

figure;
subplot(2,1,1);
stem(transmit_data);
title('transmit_data');
subplot(2,1,2);
stem(receive_data);
title('receive_data');