[信息与通信]Matlab程序文本说明.doc

Matlab 程序文本 第一章 离散时间信号与系统本程序产生均匀分布的白噪声信号set(gcf, color, w);p = 0.1; N = 500000;u = rand(1,N); % 用函数rand 产生500000 个随机数组 u_mean1 = mean(u(1:length(u) power_u = var(u) subplot(211);plot(u(1:100); grid; % 显示前100个随机数xlabel(n, fontSize, 12); ylabel(u(n), fontSize, 12); title(（1）前 100 个随机数);subplot(212)hist(u,50); grid; % 画出随机数分布的直方图xlabel(n, fontSize, 12); ylabel(u(n), fontSize, 12); title(（2）随机数分布的直方图);% u_mean = 0.4997; power_u = 0.0832 % 程序运行结果 本程序产生零均值、方差为 0.1且服从高斯分布的噪声信号set(gcf, color, w)p=0.1; N=500000;u=randn(1,N); a=sqrt(p); % 计算 a 值u=u*a; power_u=var(u); % 把 randn 得到的噪声序列乘以 a 值。subplot(221);plot(u(1:100);gridxlabel(n, fontSize, 12); ylabel(u(n), fontSize, 12); title(（1）白噪声序列);subplot(222)hist(u,50) ; grid % 画出直方图xlabel(n, fontSize, 12); ylabel(u(n), fontSize, 12); title(（2）白噪声序列的高斯分布);% power_u = 0.0997 mean_u = -8.9297e-004演示序列平移set(gcf,color,w)x = 1,2,3,4,5,4,3,2,1;m = 0:8subplot(221)H = stem(m,x); set(H, markersize, 2); gridxlabel(n, FontSize, 12); ylabel(x(n), FontSize, 12); title(（1）原序列);axis(0,15,0,6);n0=4n=m + n0; y=xsubplot(222)H = stem(n,y); set(H, markersize, 2); gridxlabel(n, FontSize, 12); ylabel(x(n + no), FontSize, 12); title(（2）平移序列);axis(0,15,0,6) 演示函数 conv 在线性卷积中的用法set(gcf, color, w);h = 4*ones(1,4),zeros(1, 4); nh=0:7;x = 1,4,3,1,zeros(1, 5); nx=-1:7;y = conv(x, h);n = -1 : 14;subplot(311); H = stem(nh, h); set(H, markersize, 6); grid; xlabel(n, FontSize, 12); ylabel(h(n), FontSize, 12); title(（1）序列 1);subplot(312); H = stem(nx, x); set(H, markersize, 6); grid; xlabel(n, FontSize, 12); ylabel(x(n), FontSize, 12); title(（2）序列 2);subplot(313); H = stem(n, y); set(H, markersize, 6); grid; xlabel(n, FontSize, 12); ylabel(y(n), FontSize, 12); title(（3）卷积所得的序列);y% y = 4 20 32 36 32 16 4 0 0 0 0 0 0 0 0 0 用Toeplitz 矩阵计算线性卷积 set(gcf,color,w);x = 1,2,3h = 3,4,5,6,7 % 给定 x 和 h 向量Nx = length(x); Nh = length(h); L=Nx+Nh-1; % 求向量长度% 生成toeplitz 矩阵 H% TOEPLITZ(C,R) is a non-symmetric Toeplitz matrix having C as its first column and R as its first row. H = toeplitz(h(1),zeros(1,Nx-1),h,zeros(1,Nx-1)y = x * H % 用向量 x 左乘 toeplitz 矩阵，求出卷积subplot(2,2,1)Hd = stem(0:L-1,y); set(Hd,markersize,5); gridxlabel(n, FontSize, 12); ylabel(y(n), FontSize, 12); title(卷积所得的序列); 演示用函数 filter 和 impz 求解差分方程set(gcf,color,w);b = 1.5; a = 1, -0.5;x = 1,zeros(1,20);y1 = filter(b, a, x)subplot(221)H=stem(0:length(y1)-1,y1); axis(-1,20,0,2); grid;xlabel(n, FontSize, 12); ylabel(y(n), FontSize, 12); title(（1）用函数 filter 求得的输出序列);set(H,markersize,3);length(y1)subplot(222)y2=impz(b,a,20);length(y2)H=stem(0:length(y2)-1,y2); axis(-1,20,0,2); grid;xlabel(n, FontSize, 12); ylabel(y(n), FontSize, 12); title(（2）用函数 impz 求得的输出序列);set(H,markersize,3); 用 FIR 和 IIR 系统分别实现五点滤波算法 set(gcf,color,w)% 产生高斯分布的噪声序列p=0.01;N=1000; % 使噪声的方差为 p=0.01u=randn(1,N); a=sqrt(p);m = mean(u)u=u*a; sigma = var(u)% 显示噪化序列（信号平均值为零）subplot(221)plot(u(1:N),k); grid; xlabel(n, FontSize, 12); ylabel(sn(k), FontSize, 12); title(（1）噪化序列, FontSize, 8) % 五点平均算法 - 使用 FIR 非因果系统subplot(222)for k = 4 : N-3 v(k) = 0.2*(u(k-2)+u(k-1)+u(k)+u(k+1)+u(k+2);endsigma_FIR1 = var(v)plot(1:N-3,v,k); grid, xlabel(n, FontSize, 12); ylabel(snFIR(n), FontSize, 12); title(（2）五点平均算法（FIR非因果系统）, FontSize, 8) % 五点平均算法 - 使用 FIR 因果系统subplot(223)for k = 5 : N-3 v(k) = 0.2*(u(k-4)+u(k-3)+u(k-2)+u(k-1)+u(k);endsigma_FIR2 = var(v)plot(1:N-3,v,k); grid, xlabel(n, FontSize, 12); ylabel(snFIR(n), FontSize, 12); title(（3）五点平均算法（FIR因果系统）, FontSize, 8) % 五点平均算法 - 使用 IIR 系统subplot(224);v(3)=0;for k = 4 : N-3 v(k) = 0.8*v(k-1)+0.2*u(k);endsigma_IIR = var(v)plot(1:N-3,v,k); grid; xlabel(n, FontSize, 12); ylabel(snIIR(n), FontSize, 12); title(（4）五点平均算法（IIR系统）, FontSize, 8); % 程序运行后在Matlab 命令窗看到的实测数据 % m =% -0.0431 % 信号的实测平均值% sigma = % 0.0089 % 信号的实测方差值 % sigma_FIR1 =% 0.0018 % 方差 ( 五点平均算法 - 使用 FIR 因果系统 )% sigma_FIR2 =% 0.0018 % 方差 ( 五点平均算法 - 使用 FIR 非因果系统 )% sigma_IIR =% 0.0010 % 方差 ( 五点平均算法 - 使用 IIR 系统 )下面计算不同算法的系统输出噪声的方差。(1) 五点平均算法 - 使用 FIR 菲因果系统系统的输出、输入噪声序列分别是和。可写为 （1）式中的、和都是随机变量。由式（1）得 （3）其中，代表平均运算。将式（3）展开后，由于假定不同随机变量不相关，故有 （4）式中右边的5项是相应随机变量的方差。但假设随机过程是平稳的，故这5项是相等的，即 （5）是输出噪声的方差，而是输入噪声的方差。由此可见，采用这种系统时，输出噪声方差 （6）（2） 五点平均算法 - 使用 FIR 因果系统 （7） 同样可以证明，采用这种系统时，输出噪声方差 （8）可见，采用五点平均算法时，不管系统是FIR因果系统或FIR非因果系统，输出噪声方差是相同的。（3） 五点平均算法 - 使用 IIR 系统 （9）故 （10）因此，采用这种系统时，输出噪声方差 (11) 以上分析结果与实测结果是吻合的。 观察离散时间信号的相关性set(gcf, color, w)p = 0.1; N = 5000; % 产生点数为 5000 的白噪声序列。其均值为零，功率为 0.1且服从高斯分布。 u=randn(1,N); a=sqrt(p);u=u*a; power_u=var(u);subplot(221); % 观察序号区间1：100的噪声序列。u1 = u(1:500);plot(u1); gridxlabel(n, FontSize, 12); ylabel(u(n), FontSize, 12); title(（1）白噪声序列);N =500u1 = u(1:N);for i =1: N % 求自相关函数（令位移为 0 49 ） u2 = u(1+i-1 : 500+i-1); a(i) = u1*u2;endsubplot(222); % 画出自相关函数k = 1:N;plot(k,a); hold on; gridaxis(-50,N,-10,60)xlabel(n, FontSize, 12); ylabel(a(n), FontSize, 12); title(（2）自相关序列); 求序列的自相关函数可以检测出序列是否含有周期成分set(gcf, color, w)p = 0.1; N = 5000; % 产生点数为 5000 的白噪声序列。其均值为零，功率为 0.1且服从高斯分布。 u=randn(1,N); a=sqrt(p);u=u*a; power_u=var(u);n =1:800s= 0.2*sin(2*pi*10/600*n+pi/2)sn = s(1:400)+u(1:400)subplot(221);plot(s(1:400), gridxlabel(n, FontSize, 12); ylabel(s(n), FontSize, 12); title(（1）正弦序列, FontSize, 8);subplot(222); plot(sn); gridxlabel(n, FontSize, 12); ylabel(sn(n), FontSize, 12); title(（2）噪化的正弦序列, FontSize, 8);subplot(223);u2 = sn(1:200);for i = 1 : 200 u3 = sn(1+i-1 : 200+i-1) a(i) = u2*u3 endk = 0:199;plot(k,a); gridxlabel(n, FontSize, 12); ylabel(a(n), FontSize, 12); title(（3）从自相关函数检出序列含有周期成分, FontSize, 8);axis(0,200,-10,10),计算矩形冲激响应序列的 DTFTset(gcf,color,w) % 置图形背景色为白N=6; nmax=32; n=0:nmax; x=ones(1,N),zeros(1,nmax+1-N); % 给出输入序列w=-9.9:0.1:9.9+1e-10;X=(sin(N*w/2)./sin(w/2).*exp(-j*(N-1)*w/2); % 进行 DTFT，给出频谱序列（关键语句）subplot(2,2,1),stem(n,x,.) % 画出输入序列axis(0,20,-0.1,1.1),grid onxlabel(n, FontSize, 12); ylabel(x(n), FontSize, 12); title(（1） 输入序列)subplot(2,2,2),plot(w,abs(X),grid on % 画出模频特性xlabel(omega (rad. / sample), FontSize, 12), ylabel(|X|, FontSize, 12); title(（2）模频特性)subplot(2,2,3),plot(w,angle(X),grid on % 画出相频特性xlabel(omega (rad. / sample), FontSize, 12); ylabel(angle of X, FontSize, 12); title(（3） 相频特性)用 Matlab 计算FT 和 DTFT 的方法。验证采样间隔越小，DTFT 就越逼近 FT。figure(1)set(gcf,color,w)% 产生连续时间信号Dt=0.00005; t=-0.005:Dt:0.005; xa=exp(-1000*abs(t);% 显示连续时间信号subplot(221); plot(t*1000,xa); grid;xlabel(t (sec), FontSize, 12); ylabel(xa(t), FontSize, 12); title(（1）连续时间信号 );% 计算傅里叶积分并显示fmax = 2000; Wmax=2*pi*fmax; K=1000; k=-K:1:K; W=k*Wmax/K; % Xa=xa*exp(-j*t*W)*Dt; % Xa 是行向量 （第12行）Xa=exp(-j*W*t)*xa*Dt; % Xa 是列向量 （第13行） subplot(222); plot(W/(2*pi),abs(Xa); grid; % 以 f(Hz) 作为频率轴xlabel(f (Hz), FontSize, 12); ylabel(|Xa(jf)|, FontSize, 12); title(（2）傅里叶积分（模值） );% 对连续时间信号采样Ts=0.0002; n=-25:1:25; x = exp(-1000*abs(n*Ts);subplot(223); H = stem(n*Ts*1000,x); set(H,markersize,2); grid; axis(-5,5,0,1.1);xlabel(n, FontSize, 12); ylabel(xa(n), FontSize, 12); title(（3）对连续时间信号采样 );% 计算 DTFT （模值）并显示它的 1 个周期fs = 1/Ts; f = 2000; K=1000; k=-K:1:K; N = 200; w=(2*pi*(f/fs)*k/K);X=x*exp(-j*n*w); % X 是行向量subplot(224); plot(w*fs/(2*pi),abs(X); grid; % 以 f(Hz) 作为频率轴 xlabel(f (Hz) , FontSize, 12); ylabel(|X(jf)| , FontSize, 12); title(（4）离散时间傅里叶变换（模值，1个周期） );% 显示计算误差% a = abs(Xa)-abs(X*Ts); % Xa 和 a 都是行向量 （第23行）a = abs(Xa)-abs(X * Ts); % a 是列向量 （第24行）error_max = max(a)% error_max = -5.2822e-006% 显示 DTFT （模值）的 4 个周期figure(2)set(gcf,color,w)fs = 1/Ts; f = 2000; K=1000; a = 5; k=-a*K:1:a*K; w=(2*pi*(f/fs)*k/K);X=x*exp(-j*n*w); subplot(211); plot(w/(2*pi),abs(X); grid; % 以归一化频率 f/fs 作为频率轴xlabel(归一化频率 , FontSize, 12); ylabel(|X(jf)|, FontSize, 12); title(（5）离散时间傅里叶变换（模值，4个周期） ); 演示相时延对信号波形的影响set(gcf, color, w)n = 0:30;fs = 1000S11=2*sin(2*pi*50*(1/fs)*n); S12=2*sin(2*pi*100*(1/fs)*n); S13=2*sin(2*pi*150*(1/fs)*n)S21=2*sin(2*pi*50*(1/fs)*n - 0.3*pi); S22=2*sin(2*pi*100*(1/fs)*n - 0.6*pi); S23=2*sin(2*pi*150*(1/fs)*n - 0.9*pi)S31=2*sin(2*pi*50*(1/fs)*n - 0.3*pi); S32=2*sin(2*pi*100*(1/fs)*n - 0.7*pi); S33=2*sin(2*pi*150*(1/fs)*n - 0.8*pi)s1 = S11 + S12 + S13; s2 = S21 + S22 + S23; s3 = S31 + S32 + S33subplot(2,2,1)plot(n,s1); grid;xlabel(n, FontSize, 12); ylabel(s1(n), FontSize, 12); title(（1）原来的合成信号 1 , fontsize, 8)subplot(2,2,2)plot(n,s2); grid;1xlabel(n, FontSize, 12); ylabel(s2(n), FontSize, 12); title(（2）合成信号 2 - 相位移与频率成正比 , fontsize, 8)subplot(2,2,3) plot(n,s3); grid; axis(1 30 -5 5);xlabel(n, FontSize, 12); ylabel(s3(n), FontSize, 12); title(（3）合成信号 3 - 相位移与频率不成正比 , fontsize, 8)演示群时延对调幅波包络线的影响set(gcf, color, w);f1 = 0.3; f2 = 0.8; f3 = 1.2; % 三个低频调制分量的相位值 fc = 10; % 载波频率 fc f1, f2, f3fs = 200; % 采样频率n = 1 : 2000;c = cos(2*pi*(fc/fs)*n) % 载波% -subplot(321);d10 = 0.3; d20 = 0.6; d30 = 1.2; % 三个低频调制分量的相位值s10 = 1 + 0.6 * cos(2*pi*(f1/fs)*n + d10); % 调幅波包络线分量 1s20 = 1 + 0.2 * cos(2*pi*(f2/fs)*n + d20); % 调幅波包络线分量 2s30 = 1 + 0.5 * cos(2*pi*(f3/fs)*n + d30); % 调幅波包络线分量 3sm0 = (s10 + s20 + s30) .* c; % 原来的调幅波plot(n, sm0); gridxlabel(n, FontSize, 12); ylabel(sm0(n), FontSize, 12); title(（1）原来的调幅波 1 )% -subplot(322);s0 = (s10 + s20 + s30); % 总的调幅波包络线plot(n, s0); gridxlabel(n, FontSize, 12); ylabel(s0(n), FontSize, 12); title(（2）原来的调幅波包络线 1 );axis(0, 2000, 2, 5);% -subplot(323);d11 = 0.9; d21 = 1.8; d31 = 3.6; % 三个低频调制分量的相位值 % 这三个值分别正比于 d10, d20, d30s11 = 1 + 0.6 * cos(2*pi*(f1/fs)*n + d11); % 调幅波包络线分量 1s21 = 1 + 0.2 * cos(2*pi*(f2/fs)*n + d21); % 调幅波包络线分量 2s31= 1 + 0.5 * cos(2*pi*(f3/fs)*n + d31); % 调幅波包络线分量 3sm1 = (s11 + s21 + s31) .* c; % 包络线无失真的调幅波plot(n, sm1); gridxlabel(n, FontSize, 12); ylabel(sm1(n), FontSize, 12); title(（3）调幅波 2 )% -subplot(324);s1 = (s11 + s21 + s31); % 总的调幅波包络线plot(n, s1); gridxlabel(n, FontSize, 12); ylabel(s1(n), FontSize, 12); title(（4）无失真的调幅波包络线 2 );axis(0, 2000, 2, 5);% -subplot(325);d12 = 0.1; d22 = -0.2; d32 = 3.5; % 三个低频调制分量的相位值 % 这三个值并不分别正比于 d10, d20, d30s12 = 1 + 0.6 * cos(2*pi*(f1/fs)*n + d12); % 调幅波包络线分量 1s22 = 1 + 0.2 * cos(2*pi*(f2/fs)*n + d22); % 调幅波包络线分量 2s32 = 1 + 0.5 * cos(2*pi*(f3/fs)*n + d32); % 调幅波包络线分量 3sm2 = (s12 + s22 + s32) .* c; % 包络线无失真的调幅波plot(n, sm2); grid xlabel(n, FontSize, 12); ylabel(sm2(n), FontSize, 12); title(（5）调幅波 3 )% -subplot(326);s2 = (s12 + s22 + s32); % 总的调幅波包络线plot(n, s2); gridxlabel(n, FontSize, 12); ylabel(s2(n), FontSize, 12); title(（6）有失真的调幅波包络线 3 );axis(0, 2000, 2, 5); 图1 程序M020804.M 的运行结果 下面介绍调幅波以及对传输系统的要求。 先给出产生调幅波的程序。运行结果示于图2.% 产生调幅波的程序%set(gcf, color, w); f1 = 0.3 % 低频调制波 fc = 10; % 载波频率 fc f1fs = 200; % 采样频率n = 1 : 2000;c = cos(2*pi*(fc/fs)*n) % 高频载波subplot(321);s0 = cos(2*pi*(f1/fs)*n);plot(n, s0); gridxlabel(n, FontSize, 14, FontWeight, Bold); ylabel(sa(n), FontSize, 14, FontWeight, Bold); title(（1）低频调制信号 , FontSize, 14, FontWeight, Bold);subplot(322);n = 1 : 500;plot(n, c(1:500); gridaxis(0, 500, -1.2, 1.2);xlabel(n, FontSize, 14, FontWeight, Bold); ylabel(sc(n), FontSize, 14, FontWeight, Bold); title(（2）高频载波 , FontSize, 14, FontWeight, Bold);subplot(323);n = 1 : 2000;s1 = 1 + 0.6 * cos(2*pi*(f1/fs)*n); % 调幅波包络线 1 sm1 = s1 .* c; % 调幅波 1n = 1 : 2000;plot(n, sm1); gridxlabel(n, FontSize, 14, FontWeight, Bold); ylabel(sm1(n), FontSize, 14, FontWeight, Bold); title(（3）调幅波： 调幅度 m 1 , FontSize, 14, FontWeight, Bold);axis(0, 2000, -2.5, 2.5); 图 2 调幅波的产生 在无线电广播和通信中，有多种方法可以有效地传送有用信息（低频），其中一种是将低频信号对高频载波进行幅度调制，使后者的幅度随前者而变。图2（1）示出低频信号。图2（2）是等幅的高频载波。设低频信号的频率和载波的频率分别为 和，而。 （1） （2）用对进行幅度调制，得调幅波为 （3）式中，为调幅度 （4）这里，是低频调制信号的幅度，是载波的幅度。图2中，已设，故，如图2（3）所示。对载波进行幅度调制后，所得的调幅波的包络线出现峰点和谷点。在调幅广播场合，高频载波的功率经过充分放大后，通过天线发射出去。在接收端，解调器对调幅波解调，取出其包络线，得到低频信号。显然，在调制时，为了使调幅波包络线不出现失真，应保证。 根据三角公式 得 （5）式（3）可以化为 （6）其中，载波（频率为）分量 （7）上旁频（频率为）分量 （8） （初相位是，是低频调制信号的初相位） 下旁频（频率为）分量 （9） （初相位是，是低频调制信号的初相位） 图3示出简单调幅波的模频特性。由图可见，简单的调幅波有3个分量，即载波和连个旁频，而旁频的幅度与调幅度有关。 图3 简单调幅波的模频特性 实际上，低频调制频率处于某个范围内，即。因此，在载波两侧形成旁频带，如图4所示。 图4 一般调幅波的模频特性 由此可知，为了使调幅波无失真地通过滤波器，这个滤波器应是带通滤波器，其模频特性在通带范围内必须是平直的，图5示出归一化的滤波器模频特性。 图5 带通滤波器的模频特性假定调幅信号频率是。则按照式（6），上、下旁频分量各有三个，即 在这种场合，为了使调幅波无失真地通通过滤波器，滤波器在通带内应具有图6所示的模频特性。而在零频到这样宽的频率范围，应具有图7所示的恒定的相时延特性。 图6 带通滤波器的模频特性 图7 带通滤波器的相时延特性 上述对相时延的要求是很难满足的。事实上也没有必要这样做。因为调幅波的包络线代表真正有用的信号，载波只不过是运载工具。在接收端，经过解调，只取出调幅波的包络线而舍去载波。因此，在考虑带通滤波器时，只需保证载波的包络线无失真即可，而无需考虑载波的波形失真。图8示出实际带通滤波器的相频特性。图中，以载波的相位为参考点，只要保证所有旁频对于该点的相位移正比于频率的变化，就能使包络线波形不会出现失真。这就是说，在滤波器通带内，相位特性的斜率必须是常数，或者说，群时延应是常数。这一点类似于对低通滤波器的要求。图9示出低通滤波器的模频特性和相频特性。为了使输出信号的波形无失真，要求模频特性是平直的，而群时延是常数。注意，这里的群时延即是相时延；但在图8中，通带中的群时延是常数，而相时延并不是常数。 图8 实际带通滤波器的相频特性 图9 低通滤波器的频率特性 演示多个余弦信号组合而成的信号的混叠效应set(gcf,color,w);t = 0:0.01:8;x = 4 + 3*cos(pi*t) + 2 * cos(2*pi*t) + cos(3*pi*t);xa = 5+5*cos(pi*t);plot(t, x, t, xa); for i = 0:12 hold on; plot(i*8/12,5+5*cos(pi*i*8/12),ro);end; xlabel(t, fontsize, 12); ylabel(x(t), fontsize, 12); title(多个余弦信号组合而成的信号的混叠效应);grid; legend(x(t),xa(t); 第二章 变换域中的离散时间系统用部分分式法进行反变换，求系统的冲激响应set(gcf, color,w);b=1; a=poly(0.9,0.9,-0.7)r,p,C=residuez(b,a) % 求极点和留数% r = 0.2461 0.5625 0.1914% p = 0.9000 0.9000 -0.7000% C = n=0:80; % 给出时间序列h=(r(1)+r(2)*p(1).n+r(2)*n