1、2022-5-301通信原理课程设计2012年12月2022-5-302通信MATLAB仿真基础信号 syms a b c solve(a*x2+b*x+c=0,x)ans = 1/2/a*(-b+(b2-4*a*c)(1/2) 1/2/a*(-b-(b2-4*a*c)(1/2)% waveform and spectrum of sin signalclose allK = 10;f1 = 1;f2 = 2;N = 2k;dt = 0.01; %msdf = 1/(N*dt); % KHzT = N*dt; % 截短时间Bs = N*df/2; % 系统带宽f = -Bs:df:Bs-df;
2、 %频域横坐标t = -T/2:dt:T/2-dt; %时域横坐标2022-5-303s=2*sin(2*pi*f1*t)+3*sin(2*pi*f2*t);f,S = T2F(t,s); % S是s的傅氏变换t,a = F2T(f,S); % a是S的傅氏反变换a=real(a);as=abs(S);f0=max(f1,f2);figuresubplot(2,1,1) %输出的频谱 plot(f,as,b);gridaxis(-2*f0,+2*f0,min(as),max(as) xlabel(f (KHz); ylabel(|S(f)| (V/KHz) subplot(2,1,2)plot
3、(t,a,black) %输出信号波形画图gridaxis(-2/f0,+2/f0,-5,5)xlabel(t(ms); ylabel(a(t)(V)2022-5-304傅里叶变换function f,sf = T2F(t,st)% calculate a signals Fourier transilationdt = t(2)-t(1);T = (t(end)-t(1)+dt);df = 1/T;N = length(st);f =(-N/2:N/2-1)*df;sf = fft(st);sf = T/N*fftshift(sf);function t,st = F2T(f,sf);% c
4、alculate the time signal using ifftdf = f(2)-f(1);Fmx = (f(end)-f(1) + df);dt = 1/Fmx;N = length(sf);T = dt*N;t = (-N/2:N/2-1)*dt;sff = fftshift(sf);st = Fmx*ifft(sff);噪声dt = 0.001; nt = 4096;t = (-nt/2:nt/2-1)*dt;n0 = 0.01;B = 2*fm;Pn = n0*B;noise = sqrt(Pn)*randn(1,length(t);2022-5-305function t s
5、t = lpf(f,sf,B)% filter an input signal using a lowpass filterdf = f(2)-f(1);T = 1/df;hf = zeros(1,length(f);bf = -floor(B/df):floor(B/df) + floor(length(f)/2);hf(bf) = 1;yf = hf.*sf;t,st = F2T(f,yf);st = real(st);低通滤波器2022-5-306%显示模拟调制的波形方法AM，文件S_AM.m%Signaldt=0.001; %时间采样间隔fmax=1; %信源最高频率fc=10; %载
6、波中心频率T=5; %信号时长N=T/dt;t=0:N-1*dt;mt=sqrt(2)*cos(2*pi*fmax*t); %信源%AM modulation A=2;s_am=(A+mt).*cos(2*pi*fc*t);%Power Spectrum Densityf,Xf=F2T (t,s_am); %调制信号频谱PSD=(abs(Xf).2)/T; %调制信号功率谱密度figure(1)subplot(211);plot(t,s_am);hold on; %画出AM信号波形plot(t,A+mt,r-); %标示AM的包络 图5-4 AM调制流图为各个变量赋初值生成信源信号m ( t
7、)对m ( t ) 进行A M 调制计算A M 信号功率谱开始画出A M 信号波形画出A M 信号功率谱结束AM调制2022-5-307title(AM调制信号及其包络); xlabel(t); subplot(212); %画出功率谱图形plot(f,PSD); axis(-2*fc 2*fc 0 1.5*max(PSD); title(AM信号功率谱);xlabel(f);2022-5-308DSB调制与解调% DSB modulation and demodulationclose allclear all% signal sourcedt = 0.001;fm = 1;fc = 10;
8、%nt = 4096;nt = 5000;t = (-nt/2:nt/2-1)*dt;T = (t(end) - t(1) + dt)%T = 5;%t = 0:dt:T;mt = sqrt(2)*cos(2*pi*fm*t); % input signal m(t), mean of m(t) is zero% DSB modulations_dsb = mt.*cos(2*pi*fc*t);% noiseN0 = 0.01;B = 2*fm;Pn = N0*B;noise = sqrt(Pn)*randn(1,length(t);s_dsb = s_dsb + noise;figure(1
9、)subplot(311)plot(t,s_dsb)hold onplot(t,mt,r)2022-5-309% DSB demodulationrt = s_dsb.*cos(2*pi*fc*t);rt = rt - mean(rt);subplot(312)plot(t,rt)hold onf,rf = T2F(t,rt);t,rt = lpf(f,rf,2*fm);subplot(313)plot(t,rt)hold onplot(t,mt/2,r)figure(2)f,sf = T2F(t,s_dsb);psf = (abs(sf).2)/T;plot(f,psf)2022-5-301
10、0 单极性非归零码是一种最简单、 最常用的基带信号形式。这种信号脉冲的零电平和正电平分别对应着二进制代码0和1，或者说，它在一个码元时间内用脉冲的有或无来对应表示0或1码。其特点是极性单一，有直流分量，脉冲之间无间隔。另外位同步信息包含在电平的转换之中，但是当出现连0或连1序列时没有位同步信息。 生成单极性非归零码的流程图如图所示 。为变量赋初值生成snrz信号画出snrz信号的波形结束开始2022-5-3011MATLAB实现程序如下： function y=snrz(x) %本函数实现将输入的一段二进制代码编为相应的单极性不归零码输出%输入x为二进制码，输出y为编好的码t0=200;t=0
11、:1/t0:length(x); %给出相应的时间序列for i=1:length(x) %计算码元的值 if x(i)=1 %如果输入信息为1 for j=1:t0 %该码元对应的点值取1 y(i-1)*t0+j)=1; end else for j=1:t0 % 如果输入信息为0，码元对应的点值取0 y(i-1)*t0+j)=0; end end end y=y,x(i); plot(t,y); 单极性非归零码单极性非归零码 为变量赋初值生成snrz信号画出snrz信号的波形结束开始2022-5-3012在命令窗口中键入x的二进制代码和函数名，就可以得到所对应的单极性不归零码输出，如输入以
12、下指令，将出现图示结果。x=1 0 1 1 0 0 1 0;snrz(x)单极性非归零码 2022-5-3013双极性非归零码的实现同单极性基本一样，只需将snrz.m中的判断得到0信息后的语句y(i-1)*t0+j)=0;中的0改为-1即可，波形图如图所示。 双极性非归零码双极性非归零码 双极性非归零码 2022-5-3014其MATLAB实现如下：（函数文件srz.m） function y=srz(x)%本函数实现将输入的一段二进制代码编为相应的单极性归零码输出%输入x为二进制码，输出y为编好的码t0=200;t=0:1/t0:length(x); %给出相应的时间序列for i=1:l
13、ength(x) %计算码元的值 if x(i)=1 %如果输入信息为1 for j=1:t0/2 y(2*i-2)*t0/2+j)=1; %定义前半段时间值为1 y(2*i-1)*t0/2+j)=0; %定义后半段时间值为0 end else for j=1:t0 %如果输入信息为0 y(i-1)*t0+j)=0; %定义所有时间值为0 end endendy=y,x(i); plot(t,y);单极性归零码单极性归零码 2022-5-3015 同理，在命令窗口中键入x的二进制代码和函数名，就可以得到所对应的单极性归零码输出，如输入以下指令，将出现图示结果。 x=1 0 1 1 0 0 1
14、0; srz(x)单极性归零码 2022-5-3016 双极性归零码的MATLAB实现同单极性也基本一样，只需将srz.m中的判断得到0信息后的语句for j=1:t0 y(i-1)*t0+j)=0;改为for j=1:t0/2 y(2*i-2)*t0/2+j)=-1; y(2*i-1)*t0/2+j)=0;即可，其波形图如图示。 双极性归零码双极性归零码 2022-5-3017设发送的数字基带信号为 ，其中 ， ， 独立同分布，+1和-1的发送概率相同，信道中加性高斯白噪声的双边功率谱密度为 ，接收滤波器为 的匹配滤波器。通过MATLAB仿真该通信系统的性能，并与理论结果对比 nsnnTtg
15、ats1, 1na 其他, 00 , 1sTttgna20N tg解：通过MATLAB采用蒙特卡罗方法仿真该通信系统误码性能的程序如下，流程图如图所示。 2022-5-3018EbN0dB=0:0.5:10;N0=10.(-EbN0dB/10);sigma=sqrt(N0/2);%理论计算的误码率Pb=0.5*erfc(sqrt(1./N0);%仿真误码率for n=1:length(EbN0dB)%产生等概信源a=sign(randn(1,100000); %离散等效接收模型 rk=a+sigma(n)*randn(1,100000); dec_a=sign(rk); %判决%计算误码率be
16、r(n)=sum(abs(a-dec_a)/2)/length(a); end为变量赋初值计算理论误码率画出理论误码率曲线结束开始计算仿真误码率画出仿真误码率曲线semilogy(EbN0dB,Pb);hold;semilogy(EbN0dB,ber,rd-);legend(理论值,仿真结果); xlabel(Eb/N0(dB);ylabel(Pb); 2022-5-30192022-5-3020% BER of digital baseband communication systemclose allSNRdB = 0:0.5:10; snr = 10.(SNRdB/10); N = 10
17、0000;% Theoretical BERBERt = 0.5*erfc(sqrt(snr);sigma = sqrt(0.5./snr);% Simulated BERfor n = 1:length(SNRdB); sig = sign(randn(1,N); noise = sigma(n)*randn(1,N); rs = sig + noise; rd = sign(rs); BERs(n) = sum(abs(sig-rd)/2)/N;endfigureplot(SNRdB,log10(BERt)hold onplot(SNRdB,log10(BERs),rs)基带系统误码率20
18、22-5-3021Ts=1;N=15;eye_num=6;a=1;N_data=1000; dt=Ts/N;t=-3*Ts:dt:3*Ts;%产生双极性数字信号d=sign(randn(1,N_data);dd=sigexpand(d,N);%基带系统冲击响应（升余弦）ht=sinc(t/Ts).*(cos(a*pi*t/Ts)./(1-4*a2*t.2/Ts2+eps);st=conv(dd,ht);tt=-3*Ts:dt:(N_data+3)*N*dt-dt;subplot(211)为变量赋初值生成双极性数字信号画出基带信号波形结束开始计算升余弦基带系统冲击响应画眼图眼图（升余弦）2022
19、-5-3022plot(tt,st);axis(0 20 -1.2 1.2);xlabel(t/Ts);ylabel(基带信号);subplot(212)%画眼图ss=zeros(1,eye_num*N);ttt=0:dt:eye_num*N*dt-dt; for k=3:50 ss=st(k*N+1:(k+eye_num)*N); drawnow; plot(ttt,ss); hold on;end;xlabel(t/Ts);ylabel(基带信号眼图); 2022-5-3023（a）接收端基带信号（b）基带信号眼图2022-5-3024眼图（理想低通）% eye of LowPasscle
20、ar allTs = 1;N_sample = 8;N = 1000;dt = Ts/N_sample;t = 0:dt:(N*N_sample-1)*dt;gt = ones(1,N_sample);d = sign(randn(1,N);a = sigexpand(d,N_sample);st = conv(a,gt);ht1 = 2.5*sinc(2.5*(t-5)/Ts);rt1 = conv(st,ht1);ht2 = sinc(t-5)/Ts);rt2 = conv(st,ht2);eyediagram(rt1,40,5)eyediagram(rt2,40,5)2022-5-302
21、5为变量赋初值生成2ASK信号开始结束画出原始二进制代码波形画出2ASK信号波形2ASK%本函数fskdigital.m实现将输入的一段二进制代码调制成相应的ask信号输出%s为输入二进制码，f为载波频率，ask为调制后输出信号t=0:2*pi/99:2*pi; m1=;c1=;for n=1:length(s)if s(n)=0; m=zeros(1,100);else s(n)=1; m=ones(1,100);endc=sin(f*t);m1=m1 m;c1=c1 cendask=c1.*m1;subplot(211);plot(m1)title(原始信号);axis(0 100*len
22、gth(s) -0.1 1.1);subplot(212);plot(ask)title(ASK信号);2022-5-30262ASK信号波形 在命令窗口中键入s的二进制代码和载波频率f，再输入函数名，就可以得到所对应的ask信号输出，其中载波频率与码元速率相同。2022-5-3027function fskdigital(s,f1,f2)%本函数fskdigital.m实现将输入的一段二进制代码调制成相应的fsk信号输出%s为输入二进制码，f1、f2分别为代码0、1对应的载波频率，fsk为调制后输出信号t=0:2*pi/99:2*pi;m1=;c1=;b1=;for n=1:length(s
23、) if s(n)=0; m=ones(1,100); c=sin(f2*t); b=zeros(1,100)2FSK2022-5-3028 else s(n)=1; m=ones(1,100); c=sin(f1*t); b=ones(1,100) end m1=m1 m; c1=c1 c; b1=b1 b;end fsk=c1.*m1; subplot(211); plot(b1,r) title(原始信号); axis(0 100*length(s) -0.1 1.1); grid on; subplot(212); plot(fsk)2022-5-30292FSK信号波形在命令窗口中键
24、入s的二进制代码和载波频率f1、f2，再输入函数名，就可以得到所对应的fsk信号输出，如输入以下指令:s=1 0 1 1 0 0 1 0; f1=200; f2=100; fskdigital 其中0信号所对应的载波频率与码元速率相同，1信号所对应的载波频率为码元速率的两倍。 2022-5-3030function pskdigital(s,f)%本函数实现将输入的一段二进制代码调制成相应的psk信号输出%s为输入二进制码，f为载波频率，psk为调制后输出信号t=0:2*pi/99:2*pi;m1=;c1=;b1=;for n=1:length(s) if s(n)=0;m=-ones(1,1
25、00); b=zeros(1,100) else s(n)=1; m=ones(1,100); b=ones(1,100) end2FSK2022-5-3031 c=sin(f*t); m1=m1 m; c1=c1 c b1=b1 b;endpsk=c1.*m1;subplot(211);plot(b1)title(原始信号);axis(0 100*length(s) -0.2 1.1);subplot(212);plot(psk);title(PSK信号);grid on;2022-5-3032在命令窗口中键入s的二进制代码和载波频率f，再输入函数名，就可以 得到所对应的psk信号输出，如输
26、入以下指令: s=1 0 1 1 0 0 1 0; f=100; pskdigital 其中载波频率与码元速率相同。 图8-17 2PSK信号波形 2022-5-3033QPSK（4PSK）性能% QPSK(4PSK)clear allfigureM = 4;snrdB = 3:0.5:10;snr = 10.(snrdB/10);E0 = 1;sigma = sqrt(1./(2*snr);error = zeros(1,length(snrdB);N = 10000;2022-5-3034for k = 1:length(snrdB) error(k) = 0; d = ceil(rand
27、(1,N)*M); s = sqrt(E0)*exp(j*2*pi/M*(d-1); r = s + sigma(k)*(randn(1,length(d) + j*randn(1,length(d); for m = 1:M rd(m,:) = abs(r-sqrt(E0)*exp(j*2*pi/M*(m-1); end for n = 1:length(s) dd1 = find(rd(:,n) = min(rd(:,n); dd(n) = dd1(1); if dd(n)=d(n) error(k) = error(k)+1; end endenderrorPe = error/N;plot(snrdB,log10(Pe),rs)hold onPet = erfc(sqrt(snr)*sin(pi/M);plot(snrdB,log10(Pet)
