通信原理实验一

实验一 信号特性分析注：I. II. III.均为实验中要调用到的m 函数一、周期信号的频谱仿真虽然Matlab 中有许多现成的频域分析工具，如fft 、ifft 等，但对通信原理的学习者来说，直接进行傅里叶变换更为直观。

为此，我们用Matlab 提供的函数为基础，编制了两个m 函数t2f.m 及t2f.m 。

t2f 是傅里叶正变换，对应 S （f ）= ⎰∞∞-)(t s dt eftj π2- (式1)f2t 是傅里叶反变换，对应 s(t)=df f S eftj ⎰∞∞-π2)( （式2）注： 式1、式2在仿真中，无限积分范围被近似在-T/2<t<T/2、-Bs<f<Bs 内。

【例一】（正弦信号的傅里叶变换）单频正弦信号的一般表达式为 s(t)=Acos(tf2π+θ)=Re{et j f)2(0θπ+}其傅里叶变换是S （f ）=A/2[)(fo f e j -δθ+)(fo f e j +-δθ]仿真程序：N=2^12; %采样点数 fs=16; %采样速度 Bs=fs/2; %系统带宽T=N/fs; %截短时间t=-T/2+[0:N-1]/fs; %时域采样点 f=-Bs+[0:N-1]/T; %频域采样点; %待观测正弦波的频率为1Hz ; %待观测正弦波的幅度为2V ; %待观测正弦波的初相为60ºs= ; %待观测的正弦波 S=t2f(s,fs); %???ss=A*real(f2t(S,fs)); %傅里叶反变换 Figure(1) %观察原始信号% 画出原始信号%在横轴标注t(ms)%在纵轴标注s(t)(V)Figure(2) %???plot(t,ss,'*') % 画出反变换后的信号 axis(0 3 -2.5 +2.5) Figure(3) % 观察幅度频谱 plot(f,abs(S)) %??? xlabel('f(khz)')ylabel('|S(f)|(V/Hz)')I ． 傅里叶正变换的m 函数 % t2f(s,fs).mfunction S=t2f(s,fs) % s 代表输入信号，S 代表s 的频谱，fs 是采样率 N=length(s); % 总样点数 T=1/fs*N; % 观察时间f=[ -N/2:(N/2-1)]/T; % 频谱采样点 tmp1=fft(s)/fs; tmp2=N*ifft(s)/fs;S(1:N/2)=tmp2(N/2+1:-1:2); S(N/2+1:N)=tmp1(1:N/2);S=S.*exp(j*pi*f*T);在这个程序中，fs 是采样率（1/△t ），N 是样点数，数组f 、s 和S 都是长为N 的数组，T 是信号的时间长度，f 是频域的采样位置，s 和S 分别是对s(t)和S （f ）的采样结果，tmp1、tmp2是中间量。

需要注意的是，直接对s 进行FFT 得到的向量对应的频域范围按归一化角频率是[0,2π),也即[0,fs)。

将右半部分[Bs,2Bs)周期性延拓到左侧的[-Bs,0)时，需要注意相位的因素。

II ． 傅里叶反变换的m 函数 % f2t(S,fs).mfunction s=f2t(S,fs)N=length(S); T=N/fs;t=[-(T/2):1/fs:(T/2-1/fs)]; % 时域采样点 tmp1=fft(S)/T;tmp2=N*ifft(S)/T;s(1:N/2)=tmp1(N/2+1:-1:2);s(N/2+1:N)=tmp2(1:N/2); s=s.*exp(-j*pi*t*fs);在这个程序中，数组t 是长为N 的数组，它是时域的采样位置。

二、功率信号的功率谱密度仿真【例二】(矩形脉冲的能量谱密度) 宽度为τ的矩形脉冲的表达式为g(t)=A0{ 22ττ<≤-t 其它其能量谱密度为Eg （f ）=|)(|2f S =]) sinc(f A [2ττ参考仿真程序：fs=64; %抽样率为fs=64khz T=64; %总观察时间为T=64ms%脉冲宽度为1ms%脉冲幅度为1VN=T*fs; dt=1/fs;t=[-T/2:dt:T/2-dt]; df=1/T;f=[-fs/2:df:fs/2-df]; g=zeros(1,N);idx=find(t>=-tau/2&t<tau/2); g(idx)=A; G=t2f(g,fs);Eg= ; % 能量谱密度Egt=(A*tau*sinc(f*tau)).^2; % 理论计算的能量谱密度 figure(1)plot(f,10*log10(Eg)) % 观察能量密度% 观测范围为-10<x<10 ,-60<y<2 figure(2)plot(f,10*log10(Egt)) %???% 观测范围为-10<x<10, -60<y<2 xlabel('kHz'); ylabel('dB/kHz')注：工程上通常习惯用分贝值观察功率谱或能量谱。

三、信号通过线性系统若线性系统的输入是x(t),输出是y （t ），则输出与输入的关系可以用卷积来描述y （t ）=⎰∞∞--τττd h t x )()( （式4），其中h （t ）是系统的单位冲激相应。

在离散时间和截短的情况下，式4对应到离散卷积hx yjjj i i∑-=仿真中更为简便的做法是借助频域关系来实现滤波Y （y ）=H(f)X(f) 【例三】（矩形脉冲通过巴特沃斯低通滤波器） 将一个宽为τ=1ms 的矩形脉冲通过一个3dB 带宽为500Hz 的6阶巴特沃思滤波器。

矩形脉冲的主瓣带宽为1khz 。

仿真中设置的时间分辨率为1/32ms ，频谱分辨率为1/64khz ，抽样率为fs=32khz ，总观察时间为T=64ms 。

参考仿真程序： fs=32;T=64;N=T*fs;dt=1/fs;t=[-T/2:dt:T/2-dt]; df=1/T;f=[-fs/2:df:fs/2-df]; g=zeros(1,N);idx=find(t>-0.5&t<=0.5);g(idx)=1;y=BTW_LPF(g,.5,fs,6); y=real(y); figure(1) plot(t,[g;y]); axis([-4,4,-0.2,1.2]); figure(2) Y=abs(t2f(y,fs)); G=abs(t2f(g,fs)); plot(f,20*log10([G ;Y]axis([-5,5,-80,0]);III. BTW_LPF(x,fm,fs,n).m 函数（巴特沃斯低通滤波器） n 阶巴特沃斯低通滤波器的传递函数为)(211)(ff Hmf nBTW+=（式5），其中fm 是滤波器的3dB 带宽。

下面是Matlab 函数将信号x （t ）通过一个3dB 带宽为fm 的巴特沃斯低通滤波器。

% BTW_LPF(x,fm,fs,n).mfunction y=BTW_LPF(x,fm,fs,n)% x 是输入，y 是输出,fm 是截止频率，fs 是采样率,n 是介数 N=length(x);T=N/fs;f=[-fs/2:1/T:fs/2-1/T]; H=(1+(f/fm)).^(2*n).^(-0.5); X=t2f(x,fs); Y=X.*H; y=f2t(Y,fs);四、白高斯噪声定义：白高斯噪声)(t n w 是持续时间无限、带宽无限的零均值平稳遍历高斯过程，其功率谱密度为 P(f2N，∞<<∞-f（式6） 仿真系统只有有效的观察时间（22T t T <≤-）,系统带宽也有限（BBssf <≤-）,因此只能仿真限带白噪声的一个时间片段n （t ）={)]([0t LPF n W22T t T <≤-其它（式7）虽然n （t ）并不是正真的白高斯噪声)(t n w ，但如果T 和Bs 足够大，可以忽略误差对所研究问题的影响。

白噪声)(t n w 的自相关函数是E[)(t n w )()](τδτ=-t nw，因此白噪声在不同时间上的采样结果形成一个独立随机序列。

仿真中对n （t ）的采样结果页必须要满足这一点。

由于对限带白噪声按奈奎斯特速率采样的结果是独立的，因此只要产生一组独立的 随机变量，就能仿真白噪声。

具体说就是N=sqrt(no*Bs)*randn(1,N);注意n （t ）的功率是BN s。

rand(1,N)产生N 个均值为0、方差为1的独立高斯随机数，再乘以NB s 0,使n 中元素的方差成为N0*Bs 。

需要特别注意的是，n 是随机过程n （t ）的一次实现的样值序列，此数组的平均功率mean （n.^2）(∑=Ni n iN121)也是一个随机变量，未必恰好等于BNs。

人们所称的“n （t ）的功率是B N s0”,是指平均功率，即E[)(2t n ]和E[∑=Ni n iN121]。

因此在仿真中，不要强制将mean （n.^2）归一化为BNs。

【例四】（白高斯噪声通过积分器） 白噪声)(t n w 的双边功率密度是Hz W N/2通过积分范围[t-Ts,t]的积分器后的输出为y(t)=dt t T t t wsn)(⎰-参考仿真程序：fs=10; % khz,时间分辨率是100us T=100; % ms ，频谱分辨率是10hz N=T*fs;dt=1/fs;t=[-T/2:dt:T/2-dt]; df=1/T;f=[-fs/2:df:fs/2-df];K=1000; % 统计平均功率谱密度时的样本数 NO=1; % 白噪声的功率谱密度，单位是W/KHz Ts=1; % 积分区间的长度是1msnTs=round(Ts/dt); % 积分区间中的样点数 n=round(Ts/dt); % 积分区间中的样点数 EP=zeros(1,N); for i=1:Knw=sqrt(NO/2*fs)*randn(1,N); % 白噪声y=zeros(1,N);for j=1:Na=max(1,j-n+1);y(j)=sum(nw(j:-1:a))*dt; % 积分endP=abs(t2f(y,fs)).^2/T;EP=EP+P;endsemilogy(f,EP/K)思考：程序中最后一句如果改成plot(f,EP/K)，图形会有什么区别？请运行对比。