2FSK数字通信系统matalb实现

2FSK数字通信系统matalb实现
最近一直想自己手动写一个通信程序，将整个通信过程一点一滴实现一下，下面将整个过程给大家分享一下下。
整个通信系统结构
英文文本作为信源，然后对信源进行ASCII编码，编码后进行了循环编码BCH（127，64），然后进行2FSK调制，再过信道，解调，信道译码、信源译码，最后准确无误得出最初的英文文本。
下面详细介绍一下整个过程。
首先是信源编码，直接对英文文本取十进制值然后将它转为二进制，每个字母用8比特表示，比较简单。
然后是BCH信道编码。编码的实现比较简单，但是这部分是我查资料最多的部分。BCH详细的编码原理这里不说了，确实比较复杂，具体实现过程为：用码字多项式乘以x^(n-k)与生成多项式相除所得余式再与他们的乘积相加。编码过程就结束了，具体实现时不一定非得在伽罗瓦域，只要模2就好了。
然后是调制部分，2FSK，也比较简单，只要对信源取反分别和载频相乘就好。
解调可以用相干解调，但在突发信号的情况下，需要极其精确的突发信号检测才能保证准确无误。解调时也可以对信号进行分帧，分帧长度取信号码元长度，做fft后检测对应频点的幅度，来解调，适用情况更多一些。这一期就现给大家分享相干解调。
信道解码部分自己实现确实比较困难，因此我调用了matalb内置函数bchdec。
信源解码比较简单，每8位转为十进制，再对应字母即可。
下面就给大家展示一下成果吧。（matalb2020b）第一个代码是2FSK调制解调，使用相干解调法。

% function Modulatation = Modulate(bin,Fs,Fc,Bw,Rb)
% 
% end
clc,clear;
Rb = 200;           %码元速率200
T = 1/Rb;           %码元宽度
Tm = 100;
M = 100;             %码元个数
totalT = T*M;       %信号持续时间

dt = T/Tm;          %采样间隔
Fs = 1/dt;
Fc = 1600;
Bw = 600;

t = 0:dt:totalT-dt;     %时间
wave = randi([0,1],1,M); %待调制数据
n = length(t);

jidai1 = zeros(1,length(t));
n = length(t);

for i = 1:length(t)
    flag = mod(i,100);
    if flag ~=0
    j = floor(i/100)+1;
    jidai1(i) = wave(j);
    else j = j+1;
    end
end

figure(1);                  % 绘制第2幅图
subplot(511);               % 窗口分割成2*1的，当前是第1个子图 
plot(t,jidai1,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形')% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,totalT,-0.1,1.1])   % 坐标范围限制

jidai2 = mod(jidai1 + 1,2);   
carry1 = cos(2*pi*(Fc+Bw/2)*t);
carry2 = cos(2*pi*(Fc-Bw/2)*t);


fsk = jidai1.*carry1 + jidai2.*carry2;

fsk=awgn(fsk,30);           % 信号fsk中加入白噪声，信噪比为SNR=20dB


fp1=2*(Fc-Bw/2-Rb);             % 第一个带通滤波器频率1
fs1=2*(Fc-Bw/2+Rb);             % 第一个带通滤波器频率2
b1=fir1(30, [fp1/Fs fs1/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数

band1=fftfilt(b1,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b1是fir滤波器的系统函数的分子多项式系数
subplot(512)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,band1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器1后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签



fp2=2*(Fc+Bw/2-Rb);             % 第二个带通滤波器频率1
fs2=2*(Fc+Bw/2+Rb);             % 第二个带通滤波器频率2
b2=fir1(30, [fp2/Fs fs2/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数

band2=fftfilt(b2,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b1是fir滤波器的系统函数的分子多项式系数
subplot(513)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,band2,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器2后的信号");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

subplot(514)  
fftband1 = abs(fft(band1))/n*2;
fftband2 = fftband1(1:n/2);
f = Fs/n-1:Fs/n:Fs/2-1;
plot(f,fftband2,'LineWidth',2) % 绘制经过低通滤波器后的信号
% axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器1后的频谱");% 标题
xlabel('频率/f');           % x轴标签
ylabel('幅度');             % y轴标签

subplot(515)  
fftband3 = abs(fft(band2))/n*2;
fftband4 = fftband3(1:n/2);
f = Fs/n-1:Fs/n:Fs/2-1;
plot(f,fftband4,'LineWidth',2) % 绘制经过低通滤波器后的信号
% axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("fsk经过带通滤波器2后的频谱");% 标题
xlabel('频率/f');           % x轴标签
ylabel('幅度');             % y轴标签



figure(2)
s1 = band1.*band1;
subplot(511)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,s1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("与低频相乘");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签


s2 = band2.*band2;
subplot(512)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,s2,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("与高频相乘");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

fp3=2*Rb;                   % 低通滤波器截止频率
b3=fir1(30, fp3/Fs, boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数

lvbo1=fftfilt(b3,s1);      % 对信号进行滤波，jt1是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
lvbo2=fftfilt(b3,s2);      % 对信号进行滤波，jt2是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
subplot(513)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,lvbo1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("低频支路");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
subplot(514)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,lvbo2,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,totalT,-1.1,1.1]);  % 设置坐标范围
title("高频支路");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签


%抽样判决
unmodulate = zeros(1,M);
for i = 1:M
    if abs(lvbo1(100*(i-1)+50))>abs(lvbo2(100*(i-1)+50))
        unmodulate(i) = 0;
    else
        unmodulate(i) = 1;
    end
end 

下面这一段是信源编码部分

 words = input("Input your string:","S") ;  %含提示的字符串输入
% words = fileattrib('"D:\兴趣\wdqwd.txt"','+w');
bits = abs(words);%将输入文本转化为ASCII码对应的10进制
data = [];
for i = 1:length(bits)
    data1 = de2bi(bits(i));
    b = data1;
    for j = 1:length(data1)
        data1(j) = b(length(data1)-j+1);
    end
        if length(data1) < 8
            data2 = [zeros(1,8-length(data1)) data1];
        else
        end
        data = [data,data2];
end

下面是信道编码部分，生成多项式gx是系统码的生成多项式，可以在许多书籍中查到，也可以是使用函数生成。

gx = 1206534025570773100045;             %其最高项次数位n-k
bin = oct2bin(gx);                       %将生成多项式用二进制数字表示
n = 127;k = 64; e = 10;
while 0 == bin(1)    
    bin = bin(2:length(bin));
end
R = mod(length(data),64);                %将待编码数据补0成为64的倍数
if R~=0
    data = [data,zeros(1,64-R)];   
end
m1 = zeros(1,64);                        %寄存分割后的信源序列
x1=zeros(1,n-k+1);               
x1(1)=1;                                 %x^n-k
encode = [];                             %最终编码数据
for i = 1:length(data)/64
    m1 = data((i-1)*64+1:(i-1)*64+64);
    c1=conv(m1,x1);              %m(x)* x^n-k    
    c1 = mod(c1,2);
    [~,r]=deconv(c1,bin);        
    r=mod(r,2);                  %取m(x)* x^n-k/g(x)的余式
    c = mod(c1+r,2);             %m(x)* x^n-k/g(x)的余式与m(x)* x^n-k相加                                     
    encode = [encode,c];
end

最后是整个完整的代码，速速Ctrl+C

% function Modulatation = Modulate(bin,Fs,Fc,Bw,Rb)
% 
% end
clc,clear;
Rb = 200;           %码元速率200
T = 1/Rb;           %码元宽度
Tm = 100;
 words = input("Input your string:","S") ;  %含提示的字符串输入
% words = fileattrib('"D:\兴趣\wdqwd.txt"','+w');
bits = abs(words);%将输入文本转化为ASCII码对应的10进制
data = [];
for i = 1:length(bits)
    data1 = de2bi(bits(i));
    b = data1;
    for j = 1:length(data1)
        data1(j) = b(length(data1)-j+1);
    end
        if length(data1) < 8
            data2 = [zeros(1,8-length(data1)) data1];
        else
        end
        data = [data,data2];
end

%%
%%信道编码BCH（127，64），纠错能力为10
gx = 1206534025570773100045;             %其最高项次数位n-k
bin = oct2bin(gx);                       %将生成多项式用二进制数字表示
n = 127;k = 64; e = 10;
while 0 == bin(1)    
    bin = bin(2:length(bin));
end
R = mod(length(data),64);                %将待编码数据补0成为64的倍数
if R~=0
    data = [data,zeros(1,64-R)];   
end
m1 = zeros(1,64);                        %寄存分割后的信源序列
x1=zeros(1,n-k+1);               
x1(1)=1;                                 %x^n-k
encode = [];                             %最终编码数据
for i = 1:length(data)/64
    m1 = data((i-1)*64+1:(i-1)*64+64);
    c1=conv(m1,x1);              %m(x)* x^n-k    
    c1 = mod(c1,2);
    [~,r]=deconv(c1,bin);        
    r=mod(r,2);                  %取m(x)* x^n-k/g(x)的余式
    c = mod(c1+r,2);             %m(x)* x^n-k/g(x)的余式与m(x)* x^n-k相加                                     
    encode = [encode,c];
end


%%
%进行2FSK调制
M = length(encode);   %信道编码后码元个数
totalT = T*M;       %信号持续时间
dt = T/Tm;          %采样间隔
Fs = 1/dt;
Fc = 1600;
Bw = 600;
t = 0:dt:totalT-dt;     %时间
n = length(t);

jidai1 = zeros(1,length(t));
n = length(t);
for i = 1:length(t)
    flag = mod(i,100);
    if flag ~=0
    j = floor(i/100)+1;
    jidai1(i) = encode(j);
    else j = j+1;
    end
end
figure(1);                  % 绘制第2幅图
subplot(511);               % 窗口分割成2*1的，当前是第1个子图 
plot(t,jidai1,'LineWidth',2) % 绘制基带码元波形，线宽为2
title('基带信号波形')% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签
axis([0,totalT,-0.1,1.1])   % 坐标范围限制

jidai2 = mod(jidai1 + 1,2);   
carry1 = cos(2*pi*(Fc+Bw/2)*t);
carry2 = cos(2*pi*(Fc-Bw/2)*t);
fsk = jidai1.*carry1 + jidai2.*carry2;

fsk=awgn(fsk,30);           % 信号fsk中加入白噪声，信噪比为SNR=20dB

fftfsk = abs(fft(fsk))/n*2;
fftfsk2 = fftfsk(1:n/2);
f = Fs/n-1:Fs/n:Fs/2-1;
subplot(512);   
plot(f,fftfsk2,'LineWidth',2) % 绘制经过低通滤波器后的信号
title("fsk频谱");% 标题
xlabel('频率/f');           % x轴标签
ylabel('幅度');             % y轴标签

%%
%解调
fp1=2*(Fc-Bw/2-Rb);             % 第一个带通滤波器频率1
fs1=2*(Fc-Bw/2+Rb);             % 第一个带通滤波器频率2
b1=fir1(30, [fp1/Fs fs1/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
band1=fftfilt(b1,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b1是fir滤波器的系统函数的分子多项式系数

fp2=2*(Fc+Bw/2-Rb);             % 第二个带通滤波器频率1
fs2=2*(Fc+Bw/2+Rb);             % 第二个带通滤波器频率2
b2=fir1(30, [fp2/Fs fs2/Fs], boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
band2=fftfilt(b2,fsk);      % 对信号进行滤波，fsk是等待滤波的信号，b1是fir滤波器的系统函数的分子多项式系数
 
%相干
s1 = band1.*band1;
s2 = band2.*band2;
%低通滤波
fp3=2*Rb;                   % 低通滤波器截止频率
b3=fir1(30, fp3/Fs, boxcar(31));% 生成fir滤波器系统函数中分子多项式的系数
lvbo1=fftfilt(b3,s1);      % 对信号进行滤波，jt1是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
lvbo2=fftfilt(b3,s2);      % 对信号进行滤波，jt2是等待滤波的信号，b3是fir滤波器的系统函数的分子多项式系数
%抽样判决
unmodulate = zeros(1,M);
for i = 1:M
    if abs(lvbo1(100*(i-1)+50))>abs(lvbo2(100*(i-1)+50))
        unmodulate(i) = 0;
    else
        unmodulate(i) = 1;
    end
end

unmodulate1 = zeros(1,length(t));
for i = 1:length(t)
    flag = mod(i,100);
    if flag ~=0
    j = floor(i/100)+1;
    unmodulate1(i) = unmodulate(j);
    else j = j+1;
    end
end
subplot(513)                % 窗口分割成3*1的，当前是第2个子图 
plot(t,unmodulate1,'LineWidth',2) % 绘制经过低通滤波器后的信号
axis([0,totalT,0,1.1]);  % 设置坐标范围
title("解调波形");% 标题
xlabel('时间/s');           % x轴标签
ylabel('幅度');             % y轴标签

%%
%信道译码
encode2 = zeros(length(unmodulate)/127,127);
for i = 1:length(unmodulate)/127
    for j = 1:127
        encode2(i,j) = unmodulate(j+127*(i-1));
    end
end

encode1 = gf(encode2);
uncode = bchdec(encode1,127,64);
undata = [];
 if mod(length(unmodulate)/127,2) == 0
    for i = 0:2:length(unmodulate)/127-1
    undata = [undata uncode(i+1,:) uncode(i+2,:)];
    end
 else
     for i = 0:2:length(unmodulate)/127-3
        undata = [undata uncode(i+1,:) uncode(i+2,:)];
        if i == length(unmodulate)-3
            break;
        end
     end
   undata = [undata uncode(length(unmodulate)/127,:)];
end
%%
%信源译码

dec = [];
for j = 1:length(undata)/8;
    c = undata(8*(j-1)+1:8*(j-1)+8);
    dec =[dec bin2dec(num2str(c.x))];
    
end
   txt = char(dec);

    fid = fopen('D:\matalb2020\uncode_text\text.txt', 'w');
    fprintf(fid, '%s', txt);
    fclose(fid);