www.pudn.com > matlab_digital_basic_transport.rar > basic.m
%信源为单极性归零码，取样点数自己定义
global dt df t f N T %定义全局变量
close all %关闭以前的应用窗口
clear Eb_N0 Pe
k=input('取样点数＝2^k, k=[13]'); %采样点数自定义，默认为2^13
if k==[], k=13; end
%Rt=input('占空比[0：1](默认0.5,为1则是NRZ码)');
%if Rt==[],Rt=0.5;
%end;
Rt=0.5; %占空比是0.5
N=2^k; %总的取样点数
L=32; %L为每个码元的取样点数
M=N/L; %M码元总数
Rb=2; %码元速率Rb=2Mb/s
Ts=1/Rb; %码元宽度Ts
dt=Ts/L; %时域取样间隔
df=1/(N*dt);%频域取样间隔
T=N*dt; %时间截断长度
Bs=N*df/2; %带宽
Na=2; %示波器扫描宽度为2个码元
alpha=0.5; %升余弦滚降系数
%f=[-Bs+df/2:df:Bs]; %频域横坐标
%t=[-T/2+dt/2:dt:T/2]; %时域横坐标
t=linspace(-T/2,T/2,N); %频域横坐标
f=linspace(-Bs,Bs,N)+eps;%时域横坐标
%。。。。。。。。。。。。。。。。。升余弦
hr1=sin(pi*t/Ts)./(pi*t/Ts);
hr2=cos(alpha*pi*t/Ts)./(1-(2*alpha*t/Ts).^2);
hr=hr1.*hr2;
HR=abs(t2f(hr));
GT=sqrt(HR);
GR=GT;
figure(1)
set(1,'Position',[10,250,500,400])
%设定窗口位置及大小
figure(2)
set(2,'Position',[515,250,500,400])
%设定窗口位置及大小
figure(3)
set(3,'Position',[360,10,300,200])
%设定窗口位置及大小
%。。。。。。。。。。。。。。。。。功率谱矩阵，初值为零
EP=zeros(size(f));
EPa=zeros(size(f));
EPr=zeros(size(f));
EPo=zeros(size(f));
for loop1=1:20
Eb_N0(loop1)=3*(loop1-1) ; %Eb/N0 in dB，横坐标
eb_n0(loop1)=10^(Eb_N0(loop1)/10);%求其功率
Eb=0.5;
n0=Eb/eb_n0(loop1); %信道的噪声谱密度
sita=n0*Bs; %信道中噪声功率
n_err=0; %误码计数
for loop2=1:10
%a=2*round(rand(1,M)); %产生随机码
a=round(rand(1,M)); %产生随机码
s=zeros(1,N);
ss=zeros(1,N);
for ii=1:Rt*Ts/dt; %Rt为占空比，dt为采样间隔，
%Rt*Ts/dt为每个码元中的非零值点数
s(ii+[0:M-1]*L)=a; %信源变成单极性归零码
end
ss(L/2:L:N)=a/dt;
SS=t2f(ss);
S=t2f(s); %信源信号的傅式变换
P=S.*conj(S)/T;
EP=(EP*(loop1-1)+P)/loop1; %信源功率谱
s1=real(f2t(SS.*GT)); %经发送滤波后的信源信号
S1=t2f(s1);
Pa=S1.*conj(S1)/T;
EPa=(EPa*(loop1-1)+Pa)/loop1;%信道前信号功率谱
n_ch=sqrt(sita)*randn(size(t)); %信道噪声
nr=real(f2t(t2f(n_ch).*GR)); %输出噪声
%dtt=Ts/2;
%C=abs(1-0.5*exp(-j*2*pi*f*dtt));%双径信道
%sr=real(f2t(S.*HR.*C))+nr; %接收信号
sr=real(f2t(SS.*HR))+nr; %接收信号
SR=t2f(sr);
Pr=SR.*conj(SR)/T;
EPr=(EPr*(loop1-1)+Pr)/ loop1;%接收信号功率谱
y=sr(L/2:L:N); %采样
aa=sign(y-0.5); %判决
%aa=sign(y-0.125); %判决
aa=(aa+1)/2; %判决输出信号
%aa=aa+1;
n_err=n_err+length(find(aa~=a)); %误码率计数
end
Pe(loop1)=n_err/(M*loop2); %误码率
e=zeros(1,N);
for jj=1:Rt*Ts/dt; %Rt为占空比，dt为采样间隔，
%Rt*Ts/dt为每个码元中的非零值点数
e(jj+[0:M-1]*L)=aa; %将判决输出信号转化为半占空比
end
E=t2f(e);
Po=E.*conj(E)/T;
EPo=(EPo*(loop1-1)+Po)/loop1;%输出信号功率谱
end
figure(1)
subplot(2,2,1)
plot(t,s,'r');
axis([-T/16,T/16,min(s)-0.5,max(s)+0.5])
xlabel('t (ms)');
ylabel('s (V/KHz)');
title('单极性归零信号波形');
subplot(2,2,2)
bb=30+10*log10(EP+eps);
plot(f,bb,'g')
axis([-20,20,-20,max(bb)+0.5])
grid
xlabel('f(Mhz)');
ylabel('Ps(f) (dBm/MHz)');
title('信源功率谱密度图形');
subplot(2,2,3)
plot(t,e,'r');
axis([-T/16,T/16,min(e)-0.5,max(e)+0.5])
xlabel('t (ms)');
ylabel('e (V/KHz)');
title('输出信号波形');
subplot(2,2,4)
bb=30+10*log10(EPo+eps);
plot(f,bb,'g')
axis([-20,20,-20,max(bb)+0.5])
grid
xlabel('f(Mhz)');
ylabel('Po(f) (dBm/MHz)');
title('输出信号功率谱密度');
figure(2)
subplot(2,3,1)
plot(t,s1,'b');
axis([-T/16,T/16,min(s1)-0.5,max(s1)+0.5])
xlabel('t (ms)');
ylabel('s1 (V/KHz)');
title('信道前信号波形');
subplot(2,3,2)
tt=[0:dt:Na*L*dt];
hold on
for jj=1:Na*L:N-Na*L
plot(tt,s1(jj:jj+Na*L)); %信道前信号眼图
end
title('信道前升余弦滚降眼图')
subplot(2,3,3)
bb=30+10*log10(EPa+eps);
plot(f,bb,'g')
axis([-20,20,-20,max(bb)+0.5])
grid
xlabel('f(Mhz)');
ylabel('Pa(f) (dBm/MHz)');
title('信道前信号功率谱密度图形');
subplot(2,3,4)
plot(t,sr,'b');
axis([-T/16,T/16,min(sr)-0.5,max(sr)+0.5])
xlabel('t (ms)');
ylabel('sr (V/KHz)');
title('接收信号波形');
subplot(2,3,5)
tt=[0:dt:Na*L*dt];
hold on
for jj=1:Na*L:N-Na*L
plot(tt,sr(jj:jj+Na*L)); %接收信号眼图
end
title('信道后升余弦滚降眼图')
subplot(2,3,6)
bb=30+10*log10(EPr+eps);
plot(f,bb,'g')
axis([-20,20,-20,max(bb)+0.5])
grid
xlabel('f(Mhz)');
ylabel('Pr(f) (dBm/MHz)');
title('接收信号功率谱密度图形');
figure(3)
semilogy(Eb_N0,Pe,'g');
eb_n0=10.^(Eb_N0/10);
hold on
plot(Eb_N0,0.5*erfc(sqrt(eb_n0/2)));
axis([0,19,1e-4,1])
xlabel('Eb/N0')
ylabel('Pe')
title('误码率曲线')