[collapse title=”仅供参考学习 抄袭与本人无关”]

沙袋豪华退学券

[/collapse]

子函数

  • 双极性信源信号产生子函数 SourceSignal

实现思路：通过rand(1,L)生成1×L，取值在(0,1)内的矩阵。当矩阵中元素的值小于0.5时，用-1替换；当大于0.5时，用1替换。即可得到取值为-1或1、长度为L的均匀分布双极性序列。

1
2
3
4
5
6
7
8
9
10
function an = SourceSignal(L)
an=rand(1,L);   %准备一个1*L，数值在(0,1)的随机序列
for i=1:L       %通过for循环判断序列的值，用1或-1替换序列的值，从而得到双极性伪随机序列
if an(i)<0.5
    an(i)=-1;   
else
    an(i)=1;
end
end
end
  • 发送信号子函数SendSignal

实现思路：首先定义一个1×A*L的矩阵dn，对矩阵dn每隔A插入an的对应元素，得到周期为A的发送信号。

1
2
3
4
5
6
7
8
%% 发送信号生成函数
function dn = SendSignal(an,A)
L = length(an); %获取序列的码元个数
dn = zeros(1,A*L);
for i = 1 : L
    dn(A*(i-1)+1)=an(i);    %插入零点
end
end
  • 升余弦滚降滤波器（非匹配）生成子函数RaisedCosine

实现思路：窗函数法设计。通过升余弦滚降滤波器的单位冲激响应时域表达式得到滤波器序列hn，但是因为在-(N-1)/2:(N-1)/2的范围内，会出现数值溢出的情况，所以需要加判断处理数溢出情况，且对对称中心赋值1。再利用布莱克曼窗的时域表达式得到窗序列wn，将hn和wn相乘完成滤波器加窗得到hf，将hf以二进制形式写入filter.bin二进制文件中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
% 抽样点数N，滚降系数a
function [hn,wn,hf] = RaisedCosine(N,a)
Tc = 4;
n = -(N-1)/2 : (N-1)/2;
% 升余弦滚降滤波器时域单位冲激响应 
hn = sin(pi*n/Tc).*cos(a*pi*n/Tc)./(pi*n/Tc)./(1-4*a*a*n.*n/Tc/Tc);
for i = 1 : N
    if(n(i) == 0)
        hn(i) = 1;
    elseif(1-4*a*a*n(i)*n(i)/Tc/Tc == 0)
        hn(i) = 0;
    end
end
% 布莱克曼窗函数表达式
n=0:N-1;
wn=0.42-0.5.*cos(2.*pi.*n./(max(n)-1))+0.08.*cos(4.*pi.*n./(max(n)-1));
% 给滤波器加窗
hf = hn .*wn;
fid = fopen('filter.bin','wb');   %以二进制数据写入方式打开文件
fwrite(fid,hf,'double');         % 将数据写入二进制文件中
fclose(fid);    %关闭文件
end
  • 平方根升余弦滤波器（匹配）生成子函数SRRaisedCosine

实现思路：频率抽样法设计。通过升余弦滚降滤波器的单位冲激响的频率响应通过if语句分段写出，得到频率响应序列Hf，再经过开方和IDFT变换得到平方根升余弦滤波器的单位冲激响应时域序列hp，将hp以二进制形式写入filter.bin二进制文件中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
%% 
function hp = SRRaisedCosine(N,a)
hp = f_sampling(N,a);
fid = fopen('filter.bin','wb');   %以二进制数据写入方式打开文件
fwrite(fid,hf,'double');         % 将数据写入二进制文件中
fclose(fid);    %关闭文件
end
%% 理想滤波器经频率抽样、IDFT后返回的平方根单位冲激响应、频响函数和频率抽样点
function hp = f_sampling(N,a)
k = -(N-1)/2:(N-1)/2;
n = k;
fs = 1;
% 抽样点确定 抽样间隔fs/N
f = k*fs/N;
Tc = 4;
f1 = (1-a)/(2*Tc);  
f2 = (1+a)/(2*Tc);
% 升余弦滚降滤波器频率抽样响应函数
Hf = zeros(1,N);
% 得到频率抽样点数值
for i = 1 : N
    if(abs(f(i))<f1)
        Hf(i) = Tc;
    elseif(abs(f(i))<f2)
        Hf(i) = Tc/2*(1+cos(pi*Tc/a*(abs(f(i))-(1-a)/(2*Tc))));
    else
        Hf(i) = 0;
    end
end
Hf1 = sqrt(Hf);
% IDFT 得到平方根升余弦数字滤波器的单位冲激响应
hp = 1/N*Hf1*exp(j*2*pi/N*k'*n);
end
  • 发送滤波器输出信号生成子函数SendOut

实现思路：将滤波器冲激响应h和发送信号dn进行卷积的xn，并对xn进行截断，最后得到发送滤波器输出信号xn。

1
2
3
4
5
6
function xn =SendOut(h,dn,A,L)
% 发送信号与发送发送滤波器单位冲激响应卷积
xn=conv(h,dn);
N=length(h);
xn=xn(round((N+1)/2):round((N+1)/2-1+L*A));
end
  • 比特能量子函数Average_energy

实现思路：在每次调用时先将Eb归0，然后遍历发送序列xn，将每个码元幅值的平方累加到Eb上，最后除以序列长度L，就得到了每比特平均能量。

1
2
3
4
5
6
7
function Eb=Average_energy(xn,L)
Eb=0;
for i = 1 : L
Eb=Eb+abs(xn(i))^2;
end
Eb=Eb/L;
end
  • 高斯噪声序列产生子函数Gaussin

实现思路：利用SNR=10lg（Eb/N0）可得出需要的噪声N0为多少，噪声的方差sigma=N0/2，利用方差*rand函数产生的标准正态分布随机数就可以得到方差为sigma的高斯分布随机数，从而实现了根据SNR产生相应大小的高斯白噪声的需求。

1
2
3
4
5
6
7
function ni=Gaussin(SNR,Eb,xn)
N0=Eb/(10^(SNR/10));
sigma=sqrt(N0/2);
ni=0+sigma*randn(1,length(xn));
histogram(ni);
title('高斯白噪声')
end
  • AWGN信号输出子函数AWMGchannel

实现思路：信道输出=信道输入+高斯白噪声

1
2
3
function aout=AWMGchannel(xn,ni)
aout=xn+ni;
end
  • 接收滤波器输出子函数receiptout

实现思路：如果使用非匹配滤波器，那么接收滤波器输出就等于信道的输出，如果使用匹配滤波器，接收滤波器的输出就等于信道输出卷积上接收滤波器的冲激响应。

1
2
3
4
5
6
7
8
function re=receiptout(N,h,L,flag,aout,A)
if flag==0
    re=aout;
else
    re=conv(aout,h);
    re=re(round((N+1)/2):round((N+1)/2-1+L*A));
end
end
  • 抽样判决子函数judgement

实现思路：眼图的实现实际上非常简单，只需要先在图上画出一个信号周期的波形（在本系统中也就是接受滤波器中的4个点），然后使用hold on，再画下一个周期的信号，如此循环直到遍历完接收滤波器输出序列就可以得到系统的眼图。

1
2
3
4
5
6
7
8
9
10
11
12
13
function [sp,jdout]=judgement(re,L,A)
jdout=zeros(L,1);
for i=1:L
   sp(i)=real(re((i-1)*A+1)); 
end
for i=1:L 
    if sp(i)>0
        jdout(i)=1;
    else
        jdout(i)=-1;
    end
end
end
  • 眼图子函数eyediagram1

实现思路：先对接收滤波器的输出进行抽样，得到长度为L的采样序列sp，然后分别判断sp的每个码元属于判决阈值的那边，从而判断发送的是哪个信号，在本系统中，那就是将sp（i）的值与0比较。

1
2
3
4
5
6
7
8
9
10
11
function eyes = eyediagram1(re,L,Tb)
eyes = 0;
figure('name','眼图');
for i = 1:(8*Tb):(Tb*(L-8))
    for j = 1:8*Tb
        eyes(j) = re(i+j-1);
    end
    plot(eyes,'b');hold on;
end
hold off;
end
  • 误比特率子函数BER

实现思路：遍历判决输出序列和信源发送序列，统计有多少码元不同，就可以得到误码率。

1
2
3
4
5
6
7
8
9
10
function BER(an,jdout,L)
BE=0;
for i=1:L
    if an(i)~=jdout(i)
    BE=BE+1;
    end
end
Pe=BE/L;
fprintf('误码率：%4.2f%%\n',Pe*100);
end
  • 星座图子函数constellation

实现思路：对于抽样序列sp，先根据发送序列是0还是1分成2堆，一堆用红色表示，另一堆用蓝色表示，然后使用plot将两堆数据花在一个图上即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function constellation(sp,L,an)
j = 1;k = 1;
   for i=1:L
     if an(i) == 1     
            x1(j)=sp(i);      
                j = j+1;
     else
            x2(k)=sp(i);       
                k = k+1;
     end
   end
figure('name','星座图'); 
plot(x1,zeros(1,length(x1)),'.b',x2,zeros(1,length(x2)),'.r');
title('星座图');
end

系统主程序

调用所有子函数实现一个数字基带传输系统的仿真，主程序的运行可视为整个数字基带系统的传输。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
close all;
clear;
clc;
% 二进制基带系统主程序
% 一个比特周期内的抽样点数 A = Tb
A = 4;
L = input("请输入信源序列长度：");
%% 信源产生
an = SourceSignal(L);   % 得到信源序列 an
dn = SendSignal(an,A);  % 得到发送信号 dn
%% 发送滤波器
flag = input("请选择滤波器模式：0. 非匹配滤波器模式   1. 匹配滤波器模式：");
a = input("请输入滤波器滚降系数α：");
N = input("请输入滤波器长度N：");
if (flag == 0)
    [hn,wn,h] = RaisedCosine(N,a);
elseif(flag == 1)
    h = SRRaisedCosine(N,a);
end
% 读入滤波器单位冲激响应数据文件
fid = fopen('filter.bin','rb');   % 以二进制数据读取方式打开文件
h = fread(fid,'double');        % 将二进制文件的数据读出
fclose(fid);
xn = SendOut(h,dn,A,L); % 得到发送滤波器输出信号
%% 求平均比特能量
Eb = Average_energy(xn,L);
%% 给定SNR生成高斯白噪声
SNR = input("请输入信噪比SNR：");
ni = Gaussin(SNR,Eb,xn);
%% 经AWGN信道
aout = AWMGchannel(xn,ni);
%% 接收滤波器
re = receiptout(N,h,L,flag,aout,A);
%% 抽样判决
[sp,jdout] = judgement(re,L,A);

滤波器性能测试与数据分析

通过仿真项目探究非匹配滤波器与匹配滤波器的性能差别和滤波器参数（滤波器长度、滚降系数）对滤波器性能（第一零点带宽、阻带最小衰减等）的影响。

滤波器性能（第一零点带宽和阻带最小衰减）测试程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
%% 测算滤波器第一零点带宽和阻带最小衰减程序
close all;
clear;
clc;
%% 测试滚降系数α
N = 40;
for a = 0.05:0.1:0.95
%% 测试滤波器长度N
% a = 0.5;
% for N = 31:2:51
%     h = SRRaisedCosine(N,a);    %测试匹配滤波器
    h = RaisedCosine(N,a);    %测试非匹配滤波器
    [Hw,w] = freqz(h);
    %% 求第一零点带宽
    wq = min(w):max(w)/length(w)/1000:max(w);
    %对幅频响应进行插值，可得到多的点，结果更为精确
    Hx = interp1(w,Hw,wq);      
    for i = 2:length(Hx)
        if((abs(Hx(i))<abs(Hx(i-1)))&&(abs(Hx(i))<abs(Hx(i+1)))&&(abs(Hx(i))<0.1))
            disp(['第一零点带宽为',num2str(wq(i)),' rad/s']);
            break;
        end
    end
    %% 求阻带最小衰减
    dbi = 20*log(abs(Hx)/max(abs(Hx)));
    for j = i:length(dbi)
        if((dbi(j)>dbi(j-1))&&(dbi(j)>dbi(j+1)))
            disp(['阻带最小衰减为',num2str(dbi(j)),' dB']);
            break;
        end
    end
    fprintf('\n');
end