信号的频域分析与系统的频域分析专题研讨汇总

信号与系统课程研究性学习手册信号的频域分析专题研讨【目的】(1) 建立工程应用中有效带宽的概念，了解有限次谐波合成信号及吉伯斯现象。(2) 掌握带限信号，带通信号、未知信号等不同特性的连续时间信号的抽样，以及抽样过程中的参数选择与确定。认识混叠误差，以及减小混叠误差的措施。(3) 加深对信号频域分析基本原理和方法的理解。(4) 锻炼学生综合利用所学理论和技术，分析与解决实际问题的能力。【研讨内容】基础题题目1：吉伯斯现象(1)以定义信号的有效带宽，试确定下图所示信号的有效带宽，取A=1，T=2。(2)画出有效带宽内有限项谐波合成的近似波形，并对结果加以讨论和比较。(3)增加谐波的项数，观察其合成的近似波形，并对结果加以讨论和比较。(a) 周期矩形信号 (b) 周期三角波信号【知识点】连续周期信号的频域分析，有效带宽，吉伯斯现象【信号频谱及有效带宽计算】（a） 周期矩形信号利用连续Fourier级数的时移特性，以及基本周期矩形脉冲信号的频谱可知该信号的频谱Cn为 取模为n=1,2,3. Cn|=0.5| Sa(0.5*n*pi)|C0|=0;|C1|=1/pi|C2|=0；|C3|=1/3*pi;P=（1*0.5*0.5+1*0.5*0.5）/2=0.25C02+2*C12+2*C22+2*C32=0.2430.225=0.9*P所以N=3；又因为该信号是奇对称的半波镜像信号，所以没有直流分量，且只有正弦分量的奇次谐波。且对于奇对称的实周期信号有an=0;bn=j*2Cn=Sa(0.5*n*pi)*sin(0.5*n*pi);（b） 周期三角波信号利用连续Fourier级数的时移特性，以及基本周期三角波信号的频谱可知该信号的频谱Cn为|C0|=0.5|C1|=2/pi*pi=0.20C02+2*C12=0.330.29=0.9*P所以N=1；又因为该信号为偶对称的实周期信号所以只有余弦分量bn=0;an=sin(0.5*n*pi)*sin(0.5*n*pi)*(-2/(n*n*pi*pi);【仿真程序】（a） 周期矩形信号%矩形波时域图及频谱图%t0=0:0.001:20; L=512;x0=1*sign(sin(pi*t0+0);figure(1);plot(t0,x0); grid on;axis(0,5,-1.5,1.5);title(时域波形);ws=10*pi;t=0:1/ws:20; x=1*sign(sin(pi*t+0);y=fft(x,L);Y=fftshift(y);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);figure(2);plot(w, abs(Y)/L); grid on;title(幅频谱图N=512);%-%矩形波基波合成%t=-20:0.001:20;N=input(N=); c0=0;XN=c0*ones(1,length(t);for n=1:2:N; XN=XN+2*sin(pi*n*t)*sinc(0.5*n)*sin(0.5*n*pi);endfigure(3);plot(t,XN); grid on;axis(-2,2,-2.2,2.2);（b）周期三角波信号%三角波时域图及频谱图%Fs=10000;t0=0:1/Fs:10;L=512;x=0.5*sawtooth(pi*t0,0.5)+0.5;figure(1);plot(t0,x); grid on;axis(0,5,-0.5,1.5);title(三角波时域波形);ws=pi;t=0:1/ws:2000; x=0.5*sawtooth(pi*t,0.5)+0.5;y=fft(x,L);Y=fftshift(y);w=(-ws/2+(0:L-1)*ws/L)/(2*pi);figure(2);plot(w, abs(Y)/L); grid on;title(三角波频谱图N=512);%-%三角波基波合成%t=-20:0.001:20;N=input(N=); c0=0.5; XN=c0*ones(1,length(t);for n=1:2:N; XN=XN+2*cos(pi*n*t)*sin(0.5*n*pi)*sin(0.5*n*pi)*(-2/(n*n*pi*pi);endfigure(3);plot(t,XN); grid on;axis(-2,2,0,1.2);【仿真结果】（a） 周期矩形信号有效带宽内N=3增加谐波次数N=10N=60（b）周期三角波信号有效带宽内N=1增加谐波次数N=10N=60【结果分析】（1）有效带宽内有限项谐波合成波形的近似度三角波要好于矩形波。（2）矩形波在时域上是奇对称的半波镜像信号，所以频域内没有直流分量，且只有正弦分量的奇次谐波。 三角波在时域上偶对称的实周期信号所以展开后余弦分量。 矩形波有跳变点，所以高频分量更多些。（3）矩形波在不连续点附近部分呈现的起伏,这个起伏的峰值大小似乎不随N 增大而下降。选取的项数越多，在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很大时，该峰起值趋于一个常数，大约等于总跳变值的9%。即产生吉伯斯现象。而三角波由于没有间断点，所以不产生吉伯斯现象，即随着N的增大，合成后的图像与原图形相似度也增大。【自主学习内容】 使用matlab生成时域信号 使用matlab进行信号的频域分析【阅读文献】1林丽莉等.信号处理与系统分析综合实验教程M.杭州 : 浙江大学出版社, 20132梁虹等.信号与线性系统分析: 基于 MATLAB 的方法与实现M.北京 : 高等教育出版社, 2006【发现问题】我们发现用理论方法推导功率较为复杂，因此，我们决定使用matlab仿真来计算有效带宽【问题探究】我们进行了以下仿真，得出了矩形的有限带宽内的谐波次数为3%用于计算功率，并确定有效带宽syms nA=1;T=2;C=;%Cn的数组C2=;%Cn平方的数组N=5;%可在此处改变谐波次数P0=(-A/2)2+(A/2)2)/T;%时域计算平均功率for n=0:N if n=0 n Cn=0; C=C Cn; C2=C2 Cn2; P=Cn2 miu=vpa(P/P0) else n Cn=A/2*sinc(n/2); C=C Cn; C2=C2,Cn2; P=P+2*Cn2 miu=vpa(P/P0)%观察效率miu的变化 endendw=0:N;stem(w,C)%画出Cn的变化情况grid on29【研讨内容】中等题题目2：分析音阶的频谱(1) 录制你所喜欢乐器（如钢琴、小提琴等）演奏的音阶，并存为wav格式。(2) 画出各音阶的时域波形，并进行比较。(3) 对所采集的音阶信号进行频谱分析，比较各音阶的频谱。【知识点】连续时间信号的频域分析【温馨提示】利用MATLAB提供的函数fft计算频谱。【题目分析】每一个音阶由基波频率不同的正弦分量构成。在非正弦的周期性振荡中，包含基波和谐波。和该振荡周期相等的正弦波分量称为基波分量。相应于这个周期的频率称为基波频率。频率等于基波频率的整倍数的正弦波分量称为谐波。由于频率不同，将产生不同的音调。【仿真程序】%读取wav文件并画声音信号的频谱y,fs=wavread(P_n5.wav);%读取P_n5.wav文件y=y(:,1);N=length(y);n=0:N-1;t=n/fs;f=n*fs/N;subplot(211);plot(t,y)xlabel(t/s);ylabel(x(t);title( P_n5音阶时域波形)Xabs=abs(fft(y,N)/N;Xabs=Xabs/max(Xabs);subplot(212);plot(f(1:N/2),Xabs(1:N/2);xlabel(f/Hz);ylabel(X(jw);title(音阶频域波形)grid on;axis(0 3000 0 1.1);%-%读取wav文件并画声音信号的频谱y,fs=wavread(P_n7.wav);%读取P_n7.wav文件y=y(:,1);N=length(y);n=0:N-1;t=n/fs;f=n*fs/N;subplot(211);plot(t,y)xlabel(t/s);ylabel(x(t);title( P_n7音阶时域波形)Xabs=abs(fft(y,N)/N;Xabs=Xabs/max(Xabs);subplot(212);plot(f(1:N/2),Xabs(1:N/2);xlabel(f/Hz);ylabel(X(jw);title(音阶频域波形)grid on;axis(0 3000 0 1.1);【仿真结果】【结果分析】提示：应从以下几方面对结果进行分析：(1) 你所选择乐器演奏的音阶，其时域波形的包络有何特点？振幅短时间快速增长后，再缓慢衰减(2) 你所选择乐器演奏的音阶，其频谱有何特点？基波是多少？谐波是多少？基波谐波do2605207801050re3005908901180mi3306509801300fa35068010401400so40079011901570la44088013001750xi49099014801980do2520105015802100【自主学习内容】1. 用MATLAB对声音信号进行时域和频域分析。2. 研究学习并比较不同音阶信号，了解了其时频波形上的区别与规律。【阅读文献】1林丽莉等.信号处理与系统分析综合实验教程M.杭州 : 浙江大学出版社, 20132梁虹等.信号与线性系统分析: 基于 MATLAB 的方法与实现M.北京 : 高等教育出版社, 2006【发现问题】(1) 改变音阶的包络，相应音阶听起来会有什么变化？(2) 音阶频谱中的谐波分量有什么作用？(3) 你所分析的乐器各音阶对应的频率是多少，之间存在什么关系？【问题探究】(1) 改变音阶的包络，相应音阶听起来会有什么变化？音色，响度变化。(2) 音阶频谱中的谐波分量有什么作用？影响着音色(3) 你所分析的乐器各音阶对应的频率是多少，之间存在什么关系？音高和频率是指数的关系其中，p是音高，f是频率标准音la,即钢琴的A4键，定义为p=69。音高每上升一个半音，p加1【研讨内容】拓展题题目3：连续时间信号的抽样(1) 对带限信号（如,等），确定合适的抽样间隔T，分析的频谱和抽样所得到离散信号的频谱X(ejW)，并将两者进行比较。 (2) 将正弦信号按抽样频率fs=8kHz进行1 秒钟抽样，得离散正弦序列xk为比较f0=2kHz, 2.2 kHz, 2.4 kHz, 2.6 kHz和 f0=7.2 kHz, 7.4 kHz, 7.6 kHz, 7.8 kHz 两组信号抽样所得离散序列的声音，解释所出现的现象。(3) 对于许多具有带通特性的信号，举例验证可否不需要满足?（1）【题目分析】抽样函数Sa(t)的频谱应该为非周期的矩形脉冲，但由于计算机处理时的信号均为离散信号，因此我们只关心其频谱在一个周期内的形状。因此，一方面，我们决定定义出Sa(t)函数后，采用Matlab的函数fft()来计算频谱。另一方面，通过信号与系统课程的学习，我们可以从理论上推导出Sa(t)函数的频谱并用计算机画图。之后将两种方法得到的频谱进行比较，观察其中的差异。最后，再增大Sa(t)函数时域的抽样间隔T，即减小抽样频率fs,观察时域与频域波形的变化。【仿真程序】%画出Sa函数时域波形及频谱的程序fx=1;fs=10;N=150*fs;%fs 取样频率;N 点数n=(0:N-1)-N/2;t=n/fs;y=sinc(fx*t/pi);%定义Sa函数subplot(311);plot(t,y);grid onxlim(-5*pi/fx 5*pi/fx);xlabel(t);ylabel(x(t);title(抽样信号Sa(t)XA=abs(fftshift(fft(y,N)/fs;%计算幅度谱subplot(312);plot(2*pi*n*fs/N,XA);xlabel(w);ylabel(X(jw);title(抽样信号Sa(t)的频谱)xlim(-2*fx 2*fx);grid on%-%理论分析得到的Sa函数的频谱temp=(pi/fx).*ones(1,fx*20);wth=linspace(-fx,fx,length(temp);subplot(313);plot(wth,temp);hold on;plot(-fx,-fx,0,pi/fx)plot(fx,fx,0,pi/fx)xlabel(w);ylabel(X(jw);title(抽样信号Sa(t)的频谱（理论）)xlim(-2*fx 2*fx);grid on%-%增大抽样间隔T后，画出Sa函数时域波形及频谱的程序fx=1;fs=0.3;N=1500;%fs 取样频率;N 点数n=(0:N-1)-N/2;t=n/fs;y=sinc(fx*t/pi);%定义Sa函数figure(2);subplot(211);plot(t,y);grid onxlim(-5*pi/fx 5*pi/fx);xlabel(t);ylabel(x(t);title(增大抽样间隔T的抽样信号Sa(t)XA=abs(fftshift(fft(y,N)/fs;%计算幅度谱subplot(212);plot(2*pi*n*fs/N,XA);xlabel(w);ylabel(X(jw);title(增大抽样间隔T的抽样信号Sa(t)的频谱)axis(-2*fx 2*fx 0 8);grid on【仿真结果】【结果分析】可以看出，当增大采样间隔T到一定值后，信号频谱将发生混叠，时域波形也发生较大失真。频域混叠后的幅度值为原来不混叠时的两倍，可知是相邻周期间的频谱发生了叠加。经过我们试验，得到当采样间隔约为1/pi时，为临界采样，及临界采样频率为pi。另外，我们发现Sa(t)函数的频域波形有非常明显的关于稳定值的震荡，可知其为吉伯斯现象。增大取样点数N后，函数在跳变点处仍然有震荡，即吉伯斯现象不消失。以下为我们画出的采样点数为15000时的频谱图，可以发现仍存在吉伯斯现象：（2）【题目分析】由Nyquist抽样定理，将正弦信号以大于其频率二倍的抽样频率抽样时，信号可以被较好还原，而当抽样频率小于正弦频率二倍抽样时，信号将发生失真。【仿真程序】%产生指定频率的正弦声音信号并按固定的抽样频率播放f0=2000 2200 2400 2600 7200 7400 7600 7800;fs=8000;N=2*fs;n=(0:N-1);t=n/fs;for i=1:8y=sin(2*pi*f0(i)*t);%定义声音信号XA=abs(fftshift(fft(y,N)/N;%计算幅度频谱f=n*fs/N;f=f-max(f)/2;sound(y,fs);subplot(211);plot(t,y);%画出时域波形xlabel(t);ylabel(x(t);grid on;title(频率f=,num2str(f0(i),Hz时域波形)xlim(0 0.02);subplot(212);plot(f,XA);%画出幅度频谱xlabel(f);ylabel(X(jw);grid on;title(频率f=,num2str(f0(i),Hz频域波形)pause(1);end【仿真结果】我们选取2000Hz和7200Hz的两个正弦信号以抽样频率fs=8000Hz抽样后的时域波形和频谱进行对比，仿真后画出的图形如下所示：【结果分析】如我们所料，7200Hz的正弦信号在以抽样频率fs=8000Hz抽样后，发生了严重的失真，由于其周期性而使得频率变为800Hz。这说明，抽样频率低于正弦信号最大频率的两倍抽样时，信号会发生失真。这种声音信号实际频率随正弦信号频率与抽样频率关系而变化的现象引起了我们的兴趣，于是我们进行了以下仿真来探索这一现象。%用来生成一个抽样频率固定，频率随时间变化的正弦声音信号fs=8000;N=0.5*fs;n=0:N-1;t=n/fs;y=sin(2*pi*0*t);for f=50:50:8000 %频率从50Hz到8000Hz间隔50Hz变化y=y sin(2*pi*f*t);endsound(y,fs)运行以上仿真程序，我们听到实际声音的频率，首先随定义的正弦信号的频率升高而升高，到4000Hz时不可闻，此时时域波形的幅度全为零。之后，实际声音的频率，随定义的正弦信号的频率升高而降低，到接近8000Hz时，音高最低。（3）【题目分析】对于许多具有带通特性的信号，其带宽远小于2fH,一般处理带通信号时常遵循带通信号采样定理，即采样频率fs2(fH-fL)即可。（其中，fH为信号的最大频率，fL为信号的最大频率）经讨论，我们决定用matlab生成一个有一定频谱特点的带限声音信号，然后对这个信号进行抽样。通过观察抽样前后的频谱变化，来验证以上分析的正确性。【仿真程序】%用来生成一个有一定频谱特点的带限声音信号fs=44100;N=0.5*fs;n=0:N-1;t=n/fs;y=sin(2*pi*1950*t);for f=2002:2:2100 %频率从2002Hz到2100Hz变化y=y sin(2*pi*f*t);endfor f=2150:2:2250 %频率从2150Hz到2250Hz变化y=y sin(2*pi*f*t);endwavwrite(y,fs,Sexy voice);%-%读取wav文件并画声音信号的频谱cleary,fs=wavread(Sexy voice.wav);y=y(:,1);N=length(y);n=0:N-1;f=n*fs/N;Xabs=abs(fft(y,N)/N;Xabs=Xabs/max(Xabs);%计算幅度频谱subplot(211);plot(f(1:N/2),Xabs(1:N/2);xlabel(f/Hz);ylabel(X(jw);title(原带限声音信号频域波形)grid on;axis(0 3000 0 1.1);fsc=646;Nc=length(y)/fs*fsc;B=N/Nc;%定义采样频率for i=1:Nc yc(i,:)=y(round(i*B),:);endNc=length(yc);nc=0:Nc-1;fc=nc*fsc/Nc;Xabsc=abs(fft(yc,Nc)/Nc;Xabsc=Xabsc/max(Xabsc);%计算幅度频谱subplot(212);plot(fc(1:Nc/2),Xabsc(1:Nc/2);xlabel(f/Hz);ylabel(X(jw);title(抽样后声音信号频域波形)grid on;axis(0 800 0 1.1);【仿真结果】【结果分析】我们生成的信号的频率从1950Hz到2250Hz之间变化，按我们学过的低通信号的Nyquist抽样定理，此时采样频率fsc必须至少大于2250*2=4500Hz，才能保证信号失真较小。但是，我们发现，当我们取采样频率为646Hz时，画出的抽样后声音信号频域波形，在形状上与原信号频域波形相差不大，这说明，我们我们可以从646Hz抽样后的信号中重构得到我们需要的原始信息中的大部分，并且失真较小。这样，便验证了带通信号采样定理的正确性。【自主学习内容】 连续信号的采样及重构 傅里叶快速算法fft的功能及使用方法 带通信号采样定理【阅读文献】1林丽莉等.信号处理与系统分析综合实验教程M.杭州 : 浙江大学出版社, 20132梁虹等.信号与线性系统分析: 基于 MATLAB 的方法与实现M.北京 : 高等教育出版社, 2006【发现问题】 若连续时间信号的最高频率未知，该如何确定对信号进行抽样的最大间隔？ 一般通过查询资料，确定自己所处理类型的信号的最大频率fm，然后，让信号通过一个截止频率为fm的低通滤波器，再依据抽样定理进行抽样。【问题探究】带通信号抽样频率确定的理论分析。带通抽样定理：一个频带限制在内的时间连续信号，信号带宽，令，这里为不大于的最大正整数。如果抽样频率满足条件， （3.1-9）则可以由抽样序列无失真的重建原始信号。系统的频域分析专题研讨【目的】(1) 加深对系统频域分析基本原理和方法的理解。(2) 加深对信号幅度调制与解调基本原理和方法的理解。(3) 锻炼学生综合利用所学理论和技术，分析与解决工程实际问题的能力。【研讨内容】题目1幅度调制和连续信号的Fourier变换本题研究莫尔斯码的幅度调制与解调。本题中信号的形式为其中信号x(t)由文件ctftmod.mat定义，可用命令Load ctftmod 将文件ctftmod.mat定义的变量装入系统内存。运行命令Load ctftmod后，装入系统的变量有af bf dash dot f1 f2 t x其中bf af： 定义了一个连续系统H(s)的分子多项式和分母多项式。可利用freqs(bf,af,w)求出该系统的频率响应，也可用sys=tf(bf,af)得到系统的模型，从而用lsim求出信号通过该系统的响应。dash dot： 给出了莫尔斯码中的基本信号dash和dot的波形f1 f2： 载波频率t: 信号x(t)的抽样点x: 信号x(t)的在抽样点上的值信号x(t)含有一段简单的消息。Agend 007的最后一句话是The future of technology lies in 还未说出最后一个字，Agend 007就昏倒了。你(Agend 008)目前的任务就是要破解Agend 007的最后一个字。该字的信息包含在信号x(t)中。信号x(t)具有式0的形式。式中的调制频率分别由变量f1和f2给出，信号m1(t)，m2(t)和m3(t)对应于字母表中的单个字母，这个字母表已用国际莫尔斯码进行编码，如下表所示：A - H O - - -V -B - I P -W -C -J - - -Q - -X -D -K -R -Y -E L -S Z -F -M - -T -G - -N -U -(1) 字母B可用莫尔斯码表示为b=dash dot dot dot，画出字母B莫尔斯码波形；(2) 用freqs(bf,af,w)画出系统的幅度响应；(3) 利用lsim求出信号dash通过由sys=tf(bf,af)定义的系统响应，解释你所获得的结果；(4)用解析法推导出下列信号的Fourier变换(5)利用(4)中的结果，设计个从x(t)中提取信号m1(t)的方案，画出m1(t)的波形并确定其所代表的字母；(6)对信号m2(t)和m3(t)重复(5)。请问Agent 008The future of technology lies in 【题目分析】我们先加载数据文件，将各变量的值读入。b为一个数组序列需要画出其波形直接用plot函数即可。第二题中freqs函数的引用本来就可以把函数的频率图画出来。第三题中要画出系统函数的响应，要知道系统函数的响应，我们要先知道dash函数的取值范围，以免出现错误。 第四题就从傅里叶变换的性质出发进行计算。第五题我们需要对信号进行解调，即频谱搬移，然后用一个低通滤波器对信号进行过滤最后保留下所需要的m(t)。通过一个一个对比字母所带表的图1,2,3信号所代表的字母。我们首先对信号x(t)进行宏观的分析，仿真程序如下：%加载数据文件，对未知信号x(t)进行初步分析clear;load ctftmodfs=length(x)/max(t);n=(0:N-1)-N/2;t=(0:N-1)/fs;subplot(211);plot(t,x);%画出时域波形xlabel(t);ylabel(x(t);grid on;title(x(t)的时域波形);XA=abs(fftshift(fft(x,N)/N;f=n*fs/N;subplot(212);plot(f,XA);%画出频谱xlabel(f);ylabel(X(jw);grid on;title(x(t)的幅度频谱);xlim(-800 800)仿真结果如下：我们可以看到x(t)的频率集中在200Hz和400Hz附近，这正是数据文件中f1、f2对应的数值，即载波频率。接下来，我们只需分析信号x(t)通过题目所给系统的响应，便可破解出最后的信息。第（4）题分析如下： 【仿真程序】（1）%画出字母b的时域波形clear;load ctftmodfs=length(x)/max(t);%确定信号x(t)的取样频率b=dash dot dot dot;N=length(b);n=0:N-1;t=n/fs;plot(t,b)xlabel(t);ylabel(b(t);grid on;title(b(t)的时域波形);（2）%画出题给系统的幅度响应和相位响应load ctftmodw=linspace(0,50*pi,100);H=freqs(bf,af,w);%系统函数subplot(211);plot(w,abs(H);grid on;xlabel(omega);ylabel(|H(jomega)| );title(系统的幅度响应);subplot(212);plot(w,angle(H);grid on;xlabel(omega);ylabel(phi(omega);title(系统的相位响应);（3）%画出dash(t)的时域波形及dash(t)通过系统后的响应load ctftmodsys=tf(bf,af);fs=length(x)/max(t);N=length(dash);n=0:N-1;t=n/fs;%fs为上面确定的x(t)的取样频率y=lsim(sys,dash,t);%求dash(t)通过系统后的响应subplot(211);plot(t,dash);grid on;xlabel(t);ylabel(dash(t);title(dash(t)原始信号)subplot(212);plot(t,y);grid on;xlabel(t);ylabel(ydash(t);title(dash(t)通过系统后的响应)（5）%解调出m1(t)后画出其原始波形load ctftmodsys=tf(bf,af);fs=length(x)/max(t);N=length(x);n=0:N-1;t=n/fs;%fs为上面确定的x(t)的取样频率m1=x.*cos(2*pi*f1*t); %解调出信号m1(t)m1o=lsim(sys,m1,t);%m1(t)通过系统plot(t,m1o);grid on;xlabel(t);ylabel(m1(t);title(m1(t)时域波形)（6）%解调出m2(t)、m3(t)后画出其原始波形load ctftmodsys=tf(bf,af);fs=length(x)/max(t);N=length(x);n=0:N-1;t=n/fs;%fs为上面确定的x(t)的取样频率m2=x.*sin(2*pi*f2*t); %解调出信号m2(