信号与系统实验指导书2014版

《信号与系统》 实

验指导书

中南大学信息科学与工程学院

2014年 3月

目 录

实验一 基本信号的生成 .......................................................................... 1 实验二 信号的基本运算 .......................................................................... 8 实验三 系统的时域分析 ........................................................................ 13 实验四 周期信号的频域分析 ................................................................ 20

1

实验一 基本信号的生成

1．实验目的

? 学会使用MATLAB产生各种常见的连续时间信号与离散时间信号； ? 通过MATLAB中的绘图工具对产生的信号进行观察，加深对常用信号的

理解；

? 熟悉MATLAB的基本操作，以及一些基本函数的使用，为以后的实验奠

定基础。

2．实验原理

⑴ 连续信号的MATLAB表示 ① 指数信号

指数信号Aeat在MATLAB中可用exp函数表示，其调用形式为

y?A?exp(a?t)

【例1】 单边衰减指数信号的MATLAB表示如下：（取A?1,a??0.4） % program exa_1.m, decaying exponential A=1;a=-0.4; t=0:0.001:10; ft=A*exp(a*t); plot(t,ft)

② 正弦信号

正弦信号Acos(?0t??)和Asin(?0t??)分别用MATLAB的内部函数cos和sin表示，其调用形式为

2

10.90.80.70.60.50.40.30.20.10012345678910A?cos(w0?t?phi)

A?sin(w0?t?phi)【例2】 正弦信号的MATLAB表示如下：（取A?1,?0?2?,???/6）

% program exa_2.m, sinusoidal signal

A=1;

w0=2*pi; phi=pi/6; t=0:0.001:8;

ft=A*sin(w0*t+phi); plot(t,ft)

10.80.60.40.20-0.2-0.4-0.6-0.8-1012345678③ 抽样函数

抽样函数Sa(t)在MATLAB中用sinc函数表示，定义为

sinc(t)?sin(?t)/(?t)

其调用形式为

y=sinc(t)

【例3】 抽样函数的MATLAB表示如下：（取t??3?3?）

% program exa_3.m, sample function

t=-3*pi:pi/100:3*pi; ft=sinc(t/pi); plot(t,ft)

10.80.60.40.20-0.2-0.4-10-8-6-4-202468103

④ 矩形脉冲信号

矩形脉冲信号在MATLAB中用rectpuls函数表示，其调用形式为

y=rectpuls(t,width)

用以产生一个幅度为1，宽度为width以t?0为对称的矩形波。width的默认值为1。

【例4】 以t?2T为对称中心的矩形脉冲信号的MATLAB表示如下：（取T?1）

% program exa_4.m, rectpuls

t=0:0.001:4; T=1;

ft=rectpuls(t-2*T,2*T); plot(t,ft)

10.90.80.70.60.50.40.30.20.1000.511.522.533.54⑤ 三角波脉冲信号

三角波脉冲信号在MATLAB中用tripuls函数表示，其调用形式为

y=tripuls(t,width,skew)

用以产生一个最大幅度为1，宽度为width的三角波。函数值的非零范围为

0时，(?width/2,width/2)；参数skew的取值范围为(?1?skew??1)；当skew?产生标准正三角波。

【例5】 三角波脉冲信号的MATLAB表示如下：

（取t??33,width?4,skew?0.5)） % program exa_5.m, tripuls

t=-3:0.001:3;

ft=tripuls(t,4,0.5); plot(t,ft)

4

10.90.80.70.60.50.40.30.20.10-3-2-10123⑵ 离散信号的MATLAB表示 ① 指数序列

指数序列Aak可用MATLAB中的数组幂运算a.?k实现。 【例6】 衰减指数序列的MATLAB表示如下：（取A?1,a??0.6） % program exa_6.m, exponential sequence k=0:10;A=1;a=-0.6; fk=A*a.^k; stem(k,fk)

10.80.60.40.20-0.2-0.4-0.6012345678910② 正弦序列

正弦序列的MATLAB表示与连续信号相同，只是用stem(k,f)画出序列的波形。 【例7】 正弦序列sin(?0k)的MATLAB实现如下：（取?0??/6）

5

% program exa_7.m, discrete-time sinusoidal signal k=0:39;

fk=sin(pi/6*k); stem(k,fk)

10.80.60.40.20-0.2-0.4-0.6-0.8-10510152025303540③ 单位脉冲序列

单位脉冲序列可借助MATLAB中的零矩阵函数zeros表示。函数zeros(1,N)产生一个由N个零组成的列向量。

【例8】 有限区间的单位脉冲序列?[k]的MATLAB实现如下：（取区间：?1010）

% program exa_8.m, delta sequence

k=-10:10;

delta=[zeros(1,10),1,zeros(1,10)]; % 或者： delta=[(k-0)==0]; stem(k,delta)

10.90.80.70.60.50.40.30.20.10-10-8-6-4-20246810④ 单位阶跃序列

单位阶跃序列可借助MATLAB中的单位矩阵函数ones表示。函数ones(1,N)产

6

生一个由N个1组成的列向量。

【例9】 有限区间的单位阶跃序列u[k]的MATLAB实现如下：（取区间：?1010）

% program exa_9.m, unit step sequence

k=-10:10;

uk=[zeros(1,10),ones(1,11)]; % 或者： uk=[k>=0]; stem(k,uk)

3．实验内容

⑴ 运行以上九个例子程序，掌握一些常用基本信号的特点及其MATLAB实现方法；改变有关参数，进一步观察信号波形的变化。 ⑵ 在 k?[?10:10] 范围内产生并画出以下信号：

a) f1[k]??[k]； b) f2[k]??[k+2]； c) f3[k]??[k-4]；

d) f4[k]?2?[k+2]??[k-4]。

⑶ 在 k?[0:31]范围内产生并画出以下信号：

a) f1[k]?sin??4k?cos??4k?； b) f2[k]?cos2??4k?； c) f3[k]?sin??4k?cos??8k?。

请问这三个信号的基波周期分别是多少？

10.90.80.70.60.50.40.30.20.10-10-8-6-4-202468107

实验二 信号的基本运算

1．实验目的

? 学会使用MATLAB完成信号的一些基本运算；

? 了解复杂信号由基本信号通过尺度变换、翻转、平移、相加、相乘、差

分、求和、微分及积分等运算来表达的方法；

? 进一步熟悉MATLAB的基本操作与编程，掌握其在信号分析中的运用特

点与使用方式。

2．实验原理

⑴ 信号的尺度变换、翻转、平移

信号的尺度变换、翻转、平移运算，实际上是函数自变量的运算。在信号的尺度变换f(at)和f[Mk]中，函数的自变量乘以一个常数，在MATLAB中可用算术运算符“﹡”来实现。在信号翻转f(?t)和f[?k]运算中，函数的自变量乘以一个负号，在MATLAB中可以直接用负号“－”写出。翻转运算在MATLAB中还可以利用fliplr(f)函数实现，而翻转后信号的坐标则可由?fliplr(k)得到。在信号时移f(t?t0)和f[k?k0]运算中，函数自变量加、减一个常数，在MATLAB中可用算术运算符“＋”或“－”来实现。

【例10】 对图示三角波f(t)，试利用MATLAB画出f(2t)和f(2?2t)的波形。

10.90.80.70.60.50.40.30.20.10-3-2-101238

解：实现f(2t)和f(2?2t)的MATLAB程序如下：

% program exa_10.m, Changed triangular pulse signal t=-3:0.001:3;

ft1=tripuls(2*t,4,0.5); subplot(2,1,1) plot(t,ft1) title('f(2t)')

ft2=tripuls((2-2*t),4,0.5); subplot(2,1,2) plot(t,ft2)

title('f(2-2t)')

程序运行结果如下图所示。

1

0.8

0.6 0.4 0.2

0-3-2

1

0.5

0 -3-2

⑵ 离散序列的差分与求和

离散序列的差分?f[k]?f[k]?f[k?1]，在MATLAB中用diff函数来实现，其调用格式为：

y?diff(f)

f(2t)-10f(2-2t)123-10123离散序列的求和?f[k]与信号相加运算不同，求和运算是把k1和k2之间的所有

k?k1k2样本f[k]加起来，在MATLAB中用sum函数来实现，其调用格式为：

9

y?sum(f(k1:k2))

【例11】 用MATLAB计算指数信号(?0.6)ku[k]的能量。 解：离散信号的能量定义为： E?lim其MATLAB程序实现如下：

% program exa_11.m, the energy of exponential sequence k=0:10;A=1;a=-0.6; fk=A*a.^k;

W=sum(abs(fk).^2) 程序运行结果为：

W =

1.5625

⑶ 连续信号的微分与积分

连续信号的微分也可以用上述diff函数来近似计算。例如，

N??k??N?Nf[k]

2y?sin(x2)??2xcos(x2)可由下列MATLAB语句来近似实现：

h=.001;x=0:h:pi;

y=diff(sin(x.^2))/h;

连续信号的定积分可由MATLAB中quad函数或quad8函数来实现。其调用格式为：

quad ('function_name',a,b)

其中function_name为被积函数名（.m文件名），a和b为指定的积分区间。

tdf(t)【例12】 用对图示三角波f(t)，试利用MATLAB画出和?f(?)d?的波形。

??dt

1

0.9

0.8

0.7

0.6

0.5

0.4 0.3 0.2 0.1

0-3-2-1012310

解：为了便于利用quad函数来计算信号的积分，将图示三角波f(t)写成MATLAB函数，函数名为functri（相应的.m文件名为functri.m），程序如下：

function yt=functri(t)

yt=tripuls(t,4,0.5);

然后利用diff和quad函数，并调用自编函数functri即可实现三角波信号f(t)的微分和积分，源程序如下：

% program exa_12_1.m, differentiation

h=0.001;t=-3:h:3;

y1=diff(functri(t))*1/h; plot(t(1:length(t)-1),y1) title('df(t)/dt') 程序运行结果如图所示：

df(t)/dt0.40.20-0.2-0.4-0.6-0.8-1-1.2-3-2-10123% program exa_12_2.m, integration t=-3:0.1:3;

for x=1:length(t)

y2(x)=quad('functri',-3,t(x)); end

plot(t,y2)

title('integral of f(t)') 程序运行结果如图所示：

integral of f(t)2.521.510.50-3-2-1012311

3．实验内容

⑴ 运行以上三个例题程序，掌握信号基本运算的MATLAB实现方法；改变有关参数，考察相应信号运算结果的变化特点与规律。 ⑵ 已知信号f(t)如下图所示：

给定信号f(t)21.81.61.41.2f(t)10.80.60.40.20-4-3-2-1 0t1234a) 用MATLAB编程复现上图； b) 画出f(2?2t)的波形；

(t)c) 画出dfdt的波形；

d) 画出?f(?)d?的波形。

??t

12

实验三 系统的时域分析

1．实验目的

? 学习并掌握连续时间系统的零状态响应、冲激响应和阶跃响应的MATLAB

求解方法；

? 学习并掌握离散时间系统的零状态响应、冲激响应和阶跃响应的MATLAB

求解方法；

? 进一步深刻理解连续时间系统和离散时间系统的系统函数零极点对系

统特性的影响；

? 学习并掌握卷积的MATLAB计算方法。 2．实验原理

⑴ 连续时间系统零状态响应的求解

LTI连续时间系统以常系数微分方程描述，系统的零状态响应可通过求解初始状态为零的微分方程得到。在MATLAB中，控制系统工具箱提供了一个用于求解零初始条件微分方程数值解的函数lsim，其调用形式为：

y=lsim(sys,f,t)

式中，t表示计算系统响应的抽样点向量，f是系统输入信号向量，sys是LTI系统模型，用来表示微分方程、差分方程、状态方程。在求解微分方程时，微分方程的LTI系统模型sys要借助tf函数获得，其调用形式为：

sys=tf(b,a)

式中，b和a分别为微分方程右端和左端各项的系数向量。例如对于三阶微分方程：

a3y???(t)?a2y??(t)?a1y?(t)?a0y(t)?b3f???(t)?b2f??(t)?b1f?(t)?b0f(t)

可用下列 MATLAB 语句：

a=[a3,a2,a1,a0]; b=[b3,b2,b1,b0]; sys=tf(b,a); 来获得LTI模型。

【例13】 如图所示力学系统中物体位移y(t)与外力f(t)的关系为：

13

d2y(t)dy(t)m?f?ksy(t)?f(t) ddt2dt

y(t) ks m f(t) fd 若该物体质量m?1kg，弹簧的弹性系数ks?100N/m，物体与地面的摩擦系数fd?2Ns/m，系统的初始储能为零，外力f(t)是振幅为10、周期为1的正弦信号，求物体的位移y(t)。

解：由已知条件，系统的输入信号为f(t)?10sin2?t，系统的微分方程为：

d2y(t)dy(t)?2?100y(t)?f(t) 2dtdt计算物体位移y(t)的MATLAB源程序如下：

% program exa_13.m, solution of differential equation

ts=0;te=5;dt=0.01;

sys=tf([1],[1 2 100]); t=ts:dt:te;

f=10*sin(2*pi*t); y=lsim(sys,f,t); plot(t,y)

xlabel('Time(sec)'); ylabel('y(t)') 程序运行结果如图所示：

y(t)-0.05-0.1-0.15-0.2-0.250.20.150.10.05000.511.5214 2.5Time(sec)33.544.55

⑵ 连续时间系统冲激响应和阶跃响应的求解

在MATLAB中, 求解系统冲激响应可应用控制系统工具箱提供的函数impulse，求解阶跃响应可利用函数step。其调用形式为：

y?impulse(sys,t)y?step(sys,t)

式中，t表示计算系统响应的抽样点向量，sys是LTI系统模型。

【例14】 在例13所述力学系统中，若外力f(t)是强度为10的冲激信号，求物体的位移y(t)。

解：由已知条件，系统的输入信号为f(t)?10?(t)，系统的微分方程为：

d2y(t)dy(t)?2?100y(t)?10?(t) dt2dt计算物体位移y(t)的MATLAB源程序如下：

% program exa_14.m, Impulse response of LTI system

ts=0;te=5;dt=0.01;

sys=tf([10],[1 2 100]); t=ts:dt:te;

y=impulse(sys,t); plot(t,y)

xlabel('Time(sec)'); ylabel('y(t)') 程序运行结果如图所示：

0.210.80.60.4

y(t)0-0.2-0.4-0.6-0.800.511.522.5Time(sec)33.544.5515

⑶ 离散时间系统零状态响应的求解

大多LTI离散时间系统都可用如下线性常系数差分方程描述：

?ay[k?i]??bii?0j?0NMjf[k?j]

其中，f[k],y[k]分别表示系统的输入和输出，N是差分方程阶数。已知差分方程的N个初始状态和输入f[k]，就可以通过编程由下式迭代计算出系统的输出：

y[k]???(ai/a0)y[k?i]??(bj/a0)f[k?j]

i?1j?0NM在零初始状态下，MATLAB信号处理工具箱提供了一个filter函数,可以计算由差分方程描述的系统的响应。其调用形式为：

y?filter(b,a,f)

式中，b?[b0,b1,b2,,bM],a?[a0,a1,a2,,aN] 分别是差分方程左、右端的系数

向量，f表示输入序列，y表示输出序列。注意，输出序列的长度与输入序列的长度相同。

【例15】 受噪声干扰的信号为f[k]?s[k]?d[k]，其中s[k]?(2k)0.9k是原始信号，d[k]是噪声。已知M点滑动平均（moving average）系统的输入输出关系为：

1y[k]?M?f[k-n]

n?0M-1试利用MATLAB编程实现用M点滑动平均系统对受噪声干扰的信号去噪。 解：系统的输入信号f[k]含有有用信号s[k]和噪声信号d[k]。噪声信号d[k]可以用rand函数产生，将其叠加在有用信号s[k]上，即得到受噪声干扰的输入信号

f[k]。对信号f[k]去噪的源程序如下（取M?5）：

% program exa_15.m, Signal Smoothing by Moving Average

Filter

R=51; % Length of input signal

% generate (-0.5,0.5) uniformly distributed random numbers

d=rand(1,R)-0.5;

16

k=0:R-1;

s=2*k.*(0.9.^k); f=s+d;

figure(1);

plot(k,d,'r-.',k,s,'b--',k,f,'g-'); xlabel('Time index k');

legend('d[k]','s[k]','f[k]'); M=5;b=ones(M,1)/M;a=1; y=filter(b,a,f); figure(2);

plot(k,s,'b--',k,y,'r-'); xlabel('Time index k'); legend('s[k]','y[k]'); 程序运行结果如图所示：

876543210-1d[k]s[k]f[k]051015202530Time index k35404550876543210s[k]y[k]051015202530Time index k3540455017

⑷ 离散时间系统单位脉冲响应的求解

在MATLAB中, 求解离散时间系统单位脉冲响应,可应用信号处理工具箱提供的函数impz，其调用形式为：

h?impz(b,a,k)

式中，b?[b0,b1,b2,,bM],a?[a0,a1,a2,,aN] 分别是差分方程左、右端的系数

向量，k表示输出序列的取值范围，h就是系统的单位脉冲响应。 【例16】 用impz函数求离散时间系统：

y[k]?3y[k?1]?2y[k?2]?f[k]

的单位脉冲响应h[k]，并与理论值 h[k]??(?1)k?2(?2)k,k?0 进行比较。 解：源程序如下：

% program exa_16.m, Impulse response of discrete system k=0:10

a=[1 3 2]; b=[1];

h=impz(b,a,k); subplot(211) stem(k,h)

title('单位脉冲响应的近似值') hk=-(-1).^k+2*(-2).^k; subplot(212) stem(k,hk)

axis(0 10 -2000 3000]) title('单位脉冲响应的理论值') 程序运行结果如图所示：

18

单位脉冲响应的近似值3000200010000-1000-2000012345678910单位脉冲响应的理论值400020000-2000012345678910⑸ 离散卷积的计算

卷积是用来计算系统零状态响应的有力工具。MATLAB信号处理工具箱提供了一个计算两个离散序列卷积和的函数conv,其调用形式为：

c?conv(a,b)

式中，a、b为待卷积两序列的向量表示，c是卷积结果。向量c的长度为向量a、

b的长度之和减1，即 length(c)=length(a)+length(b)-1。

k?【例17】 已知序列x[k]?[1,2,3,4;0,1,2,3],yk?[][1,1,1,1,1;k?,0,1,4]，

利用MATLAB编程计算x[k]?y[k]并画出卷积结果。 解：源程序如下：

% program exa_17.m, Linear Convolution

x=[1,2,3,4]; y=[1,1,1,1,1]; z=conv(x,y) N=length(z); stem(0:N-1,z) 程序运行结果：

19

z =

1 3 6 10 10 9 7 4 波形如图所示：

109876543210012345673．实验内容

⑴ 运行以上五个例题程序，掌握求解系统响应的MATLAB分析方法；改变模型参数，考察系统响应的变化特点与规律。 ⑵ 设离散系统可由下列差分方程表示：

y[k]-y[k-1]+0.9y[k-2]?f[k]

计算k?[?20:100]时的系统冲激响应。

⑶ 设h[k]?(0.9)ku(k)，输入f[k]?u[k]?u[k?10]，求系统输出y[k]?f[k]?h[k]。 （取k?[?10:50]） ⑷ 已知滤波器的传递函数：

H(z)?0.22

1?0.8z?1输入信号为f(t)?2sin(0.05?t)??(t),?(t)为随机信号。试绘出滤波器的输出信号波形。（取t?[0:100]）

20

实验四 周期信号的频域分析

1．实验目的

? 掌握周期信号傅立叶级数分解与合成的计算公式

? 掌握利用MATLAB实现周期信号傅立叶级数分解与综合方法 ? 理解并掌握周期信号频谱特点 2．实验原理

1、连续时间周期信号的分解

设有周期信号f(t)，周期为T1，角频率?1?2?f1?2?，且满足狄里赫利条件，T1则该周期信号可以展开成傅立叶级数，即可表示为一系列不同频率的正弦或复指数信号之和。傅立叶级数有三角形式和指数形式两种： 1）三角形式傅立叶级数

三角形式傅立叶级数为

f(t)?a0?a1co?s1t?b1sin?1t?a2co?s2t?b2sin?2t?..?.anco?snt?bnsin?nt?...?a0??ancosn(?1t)??bnsinn(?1t)n?1n?1?? （4-1）

式中系数按如下公式求解

1a0?T1?T12T?12f(t)dt （4-2）

2T21an???T1f(t)cosn?1tdt

T12 （4-3）

2bn?T1?T12T?12f(t)sinn?1tdt

（4-4）

2）指数形式傅立叶级数

指数形式的傅立叶级数表达式为：

f(t)??Fnen?????jn?1t,n?0,?1,?2,?3,? （4-5）

式中Fn称为傅立叶复系数，可由下式求得

21

1Fn?T1?T12T?12f(t)e?jn?1tdt （4-6）

傅立叶级数的指数形式和三角形式是等价的，其系数可相互转换，具体参见《信号与系统》教科书。

2、连续时间周期信号的傅立叶综合

任何满足狄里赫利条件的周期信号，可以表示成式（4-1）或（4-5）的和式形式，式（4-1）或（4-5）成为连续时间周期信号（CTFS）综合公式。

一般说来，傅立叶级数系数有无限个非零值，即任何具有有限个间断点的周期信号都一定有一个无限项非零系数的傅立叶级数表示。但对于数值计算来说，这是无法实现的。在实际应用中，可以用有限项傅立叶级数求和来逼近。即：

f(t)?n??N?FenNjn?1t?a0??ancos(n?1t)??bnsin(n?1t) （4-7）

n?1n?1NN当N值取得较大时，上式就是原周期信号f(t)的一个很好近似。式（4-7）常称做f(t)的截断傅立叶级数表示。

MATLAB的符号积分函数int()可以用来求解连续时间周期信号的截断傅立叶级数及傅立叶表示。

求积函数int()的具体使用格式如下：

（a） intf=int(f,v); 给出符号表达式f对指定变量v的（不带积分常数）不

定积分；

（b） intf=int(f,v,a,b); 给出符号表达式f对指定变量v的定积分； 3、利用MATLAB实现周期信号的傅立叶级数分解与综合

（1）利用MATLAB求解周期矩形脉冲傅立叶级数，并绘制出各次谐波叠加的傅立叶综合波形图。

周期矩形脉冲为f(t)??G?(t?nT)，式中??1,T?5。

n????采用三角形式傅立叶级数分解与综合形式，用式（4-2）~（4-4）求出傅立叶级数分解系数，运用MATLAB的符号运算功能，用式（4-7）实现信号的综合，谐波的阶数Nf?6。

（a）实现流程

利用MATLAB实现上述分析过程的流程如下：

? 编写子函数x=time_fun_x(t)，用符号表达式表示出周期信号在第一个周

期内的符号表达式，并赋值返回给符号变量x； ? 编写子函数y=time_fun_e(t)，求出该周期信号在绘图区间内的信号样值，

并赋值给返回变量y；

? 编写求解信号傅立叶系数及绘制合成波形图的通用CTFShchsym.m，该

函数流程如下：

① 调用函数time_fun_x(t)，获取周期信号的符号表达式； ② 求出信号的傅立叶系数； ③ 求出各次谐波；

④ 绘制各次谐波叠加波形图；

22

⑤ 调用函数time_fun_e(t)，绘制原信号波形图。

（b）MATLAB算法提示及说明

? 采用符号积分int求一个周期内时间函数的三角级数展开系数：a0?A0，

an?As，bn?Bs，即计算式（4-2）~（4-4）的值；

? 用循环语句for…end求出三角级数展开系数an，bn的数值，分别为A_sym，B_sym；

? 用disp()语句输出三角级数展开系数A_sym，B_sym； ? 用傅立叶三角级数展开式（4-7）合成连续时间信号；

? 化简表达式，据函数奇偶性可知，若f(t)为奇函数，则an?0；若f(t)为

偶函数，则bn?0。

（c）源程序

编写函数文件CTFShchsym.m，这是一个计算连续时间周期信号的三角级数前6次展开系数，再用这6次谐波合成原连续时间周期信号的程序，如下所示。 function [A_sym,B_sym]=CTFShchsym

% 采用符号计算求一个周期内连续时间函数f的三角级数展开系数,再用这些 % 展开系数合成连续时间函数f.傅立叶级数 % 函数的输入输出都是数值量 % Nf=6 谐波的阶数 % Nn 输出数据的准确位数

% A_sym 第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数 % B_sym 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 % tao=1 tao/T=0.2 syms t n k x T=5;tao=0.2*T;a=0.5; if nargin<4;Nf=6;end if nargin<5;Nn=32;end x=time_fun_x(t);

A0=int(x,t,-a,T-a)/T; %求出三角函数展开系数A0

As=int(2*x*cos(2*pi*n*t/T)/T,t,-a,T-a); %求出三角函数展开系数As Bs=int(2*x*sin(2*pi*n*t/T)/T,t,-a,T-a); %求出三角函数展开系数Bs

A_sym(1)=double(vpa(A0,Nn)); %获取串数组A0所对应的ASC2码数值数组 for k=1:Nf

A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); %获取串数组A所对应的ASC2码数值数组

B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); %获取串数组B所对应的ASC2码数值数组 end

23

if nargout==0

c=A_sym;disp(c) %输出c为三角级数展开系数:第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数

d=B_sym;disp(d) %输出d为三角级数展开系数: 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 t=-8*a:0.01:T-a;

f1=c(1)+c(2).*cos(2*pi*1*t/5)+0.*sin(2*pi*1*t/5); ; % 基波

f2=c(3).*cos(2*pi*2*t/5)+0.*sin(2*pi*2*t/5); ; % 2次谐波 f3=c(4).*cos(2*pi*3*t/5)+0.*sin(2*pi*3*t/5); % 3次谐波 f4=c(5).*cos(2*pi*4*t/5)+0.*sin(2*pi*4*t/5); ; % 4次谐波 f5=c(6).*cos(2*pi*6*t/5)+0.*sin(2*pi*6*t/5); % 6次谐波 f6=f1+f2; % 基波+2次谐波

f7=f6+f3; % 基波+2次谐波+3次谐波

f8=f7+f4+f5; % 基波+2次谐波+3次谐波+4次谐波+6次谐波 subplot(2,2,1) plot(t,f1),hold on

y=time_fun_e(t) %调用连续时间函数-周期矩形脉冲 plot(t,y,'r:')

title('周期矩形波的形成—基波+直流') axis([-4,4.5,-0.5,1.5]) subplot(2,2,2) plot(t,f6),hold on y=time_fun_e(t) plot(t,y,'r:')

title('周期矩形波的形成—基波+2次谐波') axis([-4,4.5,-0.5,1.5]) subplot(2,2,3) plot(t,f7),hold on y=time_fun_e(t) plot(t,y,'r:')

title('基波+2次谐波+3次谐波') axis([-4,4.5,-0.5,1.5]) subplot(2,2,4) plot(t,f8),hold on y=time_fun_e(t) plot(t,y,'r:')

title('基波+2次谐波+3次谐波+4次谐波+6次谐波') axis([-4,4.5,-0.5,1.5]) end

%------------------------------------------- function x=time_fun_x(t)

% 该函数是CTFShchsym.m的子函数。它由符号变量和表达式写成。 h=1;

24

x1=sym('Heaviside(t+0.5)')*h; x=x1-sym('Heaviside(t-0.5)')*h;

%------------------------------------------- function y=time_fun_e(t)

% 该函数是CTFShchsym.m的子函它由符号函数和表达式写成 a=0.5;T=5;h=1;tao=0.2*T;t=-8*a:0.01:T-a; e1=1/2+1/2.*sign(t+tao/2); e2=1/2+1/2.*sign(t-tao/2);

y=h.*(e1-e2); %连续时间函数-周期矩形脉冲 （d）程序运行结果

在MATLAB命令窗口键入CTFShchsym，并按回车键，即可绘制出周期矩形波形各次谐波合成波形，如图4-1所示。

图4-1 周期矩形脉冲傅立叶分解与综合结果

执行如下命令就可以得到三角级数展开系数： [A_sym,B_sym]=CTFShchsym 结果为 A_sym =

0.2000 0.3742 0.3027 0.2018 0.0935 0.0000 -0.0624

B_sym =

0 0 0 0 0 0 0 4、周期信号频谱分析

对周期为T1 的信号f(t)进行傅立叶级数展开可得到

f(t)??Fnen????jn?1t?a0??ancos(n?1t)??bnsin(n?1t)

n?1n?1??其中

25

??12? T11(an?jbn) （3-8） 2Fn?Fnej?n?如果求出an和bn，根据以下两式可以画出周期信号的幅度谱Fn~?和相位谱?n~?。

Fn?12an?bn2 幅频 2???arctan 相频

nbnan5、MATLAB实现

由上述分析可知，只要求出了周期信号的傅立叶级数幅度Fn和相位?n，就可以根据Fn和?n随频率??n?1的变化关系画出幅度谱和相位谱。

现仍以周期矩形脉冲为f(t)??G?(t?nT)，（??1,T?5）为例，来说明如

n????何用MATLAB绘制周期信号的频谱图，并对周期信号的频谱特性进行分析。 由于绘制频谱图的前提是必须先求出周期信号的傅立叶级数系数，因此只需对上述给出的求周期信号傅立叶级数的函数CTFShchsym.m进行适当修改，即可编写出绘制周期信号频谱的通用函数。

需注意的是，由于周期信号的频谱是离散的，故在绘制频谱时，采用的是stem命令而不是plot命令，下面给出实现上述通用程序CTFStpshsym.m的过程。 （1）实现流程

此处采用三角形式傅立叶级数分解，求出分解系数an和bn，再利用式（3-8）求出傅立叶复指数系数Fn，画出Fn的幅度谱Fn~?。谐波的阶数Nf可任意指定，此处指定Nf?60。

（a）实现流程如下：

? 编写子函数y=time_fun_s(t)，用符号表达式表示出周期信号在第一个周

期内的符号表达式，并赋值返回给符号变量y； ? 编写子函数x=time_fun_e(t)，求出该周期信号在绘图区间内的信号样值，

并赋值给返回变量x；

? 编写求解信号傅立叶复指数系数Fn的通用函数，该函数流程如下：

① 调用函数time_fun_s(t)，获取周期信号的符号表达式； ② 求出信号的三角形式的傅立叶系数an和bn；

26

③ 求出信号的傅立复指数系数Fn； ④ 绘制Fn的幅度频谱图；

⑤ 调用函数time_fun_e(t)绘制信号波形。

（b）算法提示

? 采用符号积分int求[0,T]内时间信号的三角级数展开系数：a0?A0，

an?As，bn?Bs;

? 用循环语句for…end求出三角级数展开系数an和bn的数值，并赋值给变

量A_sym(k+1)， B_sym(k+1)；

? 从三角级数展开系数an和bn得到复指数展开系数Fn：（3-8）Fn可以根据式

求出，需要注意的是an、bn和Fn的自变量取值情况，an、bn的变量n的取值范围为n?0,1,2,3,?,N，而Fn的变量n的取值范围为

n?0,?1,?2,?3,?,?N，为了从an和bn得

到Fn，需要用到MATLAB的反折函数fliplr来实现频谱的反折；例如已知一个

序列a=[5 4 3 2 1]，采样位置为：0，1，2，3，4。要组成员序列反折后再叠加原序列的一个新序列（采样位置：-4，-3，-2，-1，0，1，2，3，4）。程序如下：n=0:4;

a=[5 4 3 2 1];

subplot(2,1,1);stem(n,a); b=fliplr(a); k=-4:4;

c=[b,a(2:end)];

subplot(2,1,2);stem(k,c);

结果如右图所示。

（c）源程序代码

编写函数文件CTFStpshsym.m，如下所示： [CTFStpshsym.m]

function [A_sym,B_sym]=CTFSshbpsym(T,Nf)

% 采用符号计算求[0,T]内时间函数的三角级数展开系数。 % 函数的输入输出都是数值量 % Nn 输出数据的准确位数

27

% A_sym 第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数

% B_sym 第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数 % T T=m*tao， 信号周期 % Nf 谐波的阶数 % Nn 输出数据的准确位数

% m (m=T/tao)周期与脉冲宽度之比,如m=4,8,16,100等 % tao 脉宽:tao=T/m syms t n y

if nargin<3;Nf=input('please Input 所需展开的最高谐波次数:Nf=');end T=input('please Input 信号的周期T='); if nargin<5;Nn=32;end y=time_fun_s(t); A0=2*int(y,t,0,T)/T;

As=int(2*y*cos(2*pi*n*t/T)/T,t,0,T); Bs=int(2*y*sin(2*pi*n*t/T)/T,t,0,T); A_sym(1)=double(vpa(A0,Nn)); for k=1:Nf

A_sym(k+1)=double(vpa(subs(As,n,k),Nn)); B_sym(k+1)=double(vpa(subs(Bs,n,k),Nn)); end

if nargout==0

S1=fliplr(A_sym); %对A_sym阵左右对称交换

S1(1,k+1)=A_sym(1); %A_sym的1*k阵扩展为1*(k+1)阵

S2=fliplr(1/2*S1); %对扩展后的S1阵左右对称交换回原位置 S3=fliplr(1/2*B_sym); %对B_sym阵左右对称交换 S3(1,k+1)=0; %B_sym的1*k阵扩展为1*(k+1)阵

S4=fliplr(S3); %对扩展后的S3阵左右对称交换回原位置

S5=S2-i*S4; % 用三角函数展开系数A、B值合成付里叶指数系数 S6=fliplr(S5); N=Nf*2*pi/T; k2=-N:2*pi/T:N; S7=[S6,S5(2:end)]; subplot(3,3,3)

x=time_fun_e(t) % 调用连续时间函数-周期矩形脉冲 subplot(2,1,1)

stem(k2,abs(S7)); %画出周期矩形脉冲的频谱（T=M*tao） title('连续时间函数周期矩形脉冲的双边幅度谱') axis([-80,80,0,0.12]) line([-80,80],[0,0]) line([0,0],[0,0.12]) end

28

%------------------------------------------- function y=time_fun_s(t)

% 该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。 syms a a1

T=input('please Input 信号的周期T='); M=input('周期与脉冲宽度之比M='); A=1;tao=T/M;a=tao/2;

y1=sym('Heaviside(t+a1)')*A; y=y1-sym('Heaviside(t-a1)')*A; y=subs(y,a1,a); y=simple(y);

%------------------------------ function x=time_fun_e(t)

% 该函数是CTFSshbpsym.m的子函数。它由符号变量和表达式写成。 % t 是时间数组 % T 是周期 duty=tao/T=0.2 T=5;t=-2*T:0.01:2*T;tao=T/5;

x=rectpuls(t,1); %产生一个宽度tao=1的矩形脉冲 subplot(2,2,2) plot(t,x) hold on

x=rectpuls(t-5,1); %产生一个宽度tao=1的矩形脉,中心位置在t=5处 plot(t,x) hold on

x=rectpuls(t+5,1); %产生一个宽度tao=1的矩形脉,中心位置在t=-5处 plot(t,x)

title('周期为T=5，脉宽tao=1的矩形脉冲') axis([-10,10,0,1.2]) （d）运行结果及分析 ? 程序运行结果

调用CTFStpshsym.m函数文件，即可绘出周期矩形脉冲波形信号的双边幅度频谱。指令如下：在MATLAB命令窗口键入CTFStpshsym，回车，命令窗口将会出现：please Input所需展开的最高次谐波次数Nf=?；此处输入Nf=60，然后命令窗口将出现：please Input信号的周期T=?，此处输入T=5；然后命令窗口又出现：周期与脉冲宽度之比M=?，此处输入M=5，即可绘出频谱图，如图4-2所示。

29 图3-2 周期信号的双边频谱图

由图4-2可见，周期信号的频谱是离散谱，仅含有??n?1的各分量，其相邻谱线的间隔是?1。

? 脉冲宽度与频谱的关系

脉冲与频谱关系如图4-3所示。

（a）周期T=5，脉冲宽度tao=T/4=1.25的矩形脉冲频谱

（b）周期T=5，脉冲宽度tao=T/8=0.625的矩形脉冲频谱

（c）周期T=5，脉冲宽度tao=T/16=0.3125的矩形脉冲频谱

图4-3 脉冲宽度与频谱关系

30

由图4-3可见，由于周期T相同，因而相邻谱线的间隔相同；脉冲宽度?越窄，其频谱包络线第一个零点的频率越高，即信号脉宽越宽，这就验证了信号的频带宽度与脉冲宽度?成反比，即

2?B??

?

? 周期与频谱关系

周期与频谱的关系如图4-4所示。

（a）脉冲宽度tao=1，周期T=4*tao=4的矩形脉冲频谱

（b）脉冲宽度tao=1，周期T=8*tao=8的矩形脉冲频谱

（c）脉冲宽度tao=1，周期T=16*tao=16的矩形脉冲频谱

图4-4 周期与频谱的关系

由图3-4可见，由于周期信号的时域脉冲宽度不变，这时频谱包络线的零点所在位置不变，而当周期变长时，相邻谱线的间隔减少，频谱变密。如果周期无限长（这时信号变为非周期），那么，相邻谱线的间隔将趋近于零，周期信号的频谱就过渡到非周期信号的连续谱。随着周期信号的增长，各谐波分量的幅度也相应减小。

31

3. 实验内容

1.仿照例程，实现下述周期信号的傅立叶级数分解与合成：

f(t)

1

O 4 -4 -3 1

要求：

5 t （a）首先，推导出求解a0，an，bn的公式，计算出前10次系数； （b）利用MATLAB求解a0，an，bn的值，其中an，bn求解前10次系

数，并给出利用这些系数合成的信号波形。

syms t n;

ft=0;T=4;w=2*pi/T; a0=int(2/T,t,0,1);

an=int(2*cos(n*t*w)/T,t,0,1)+int(0*2*cos(n*t*w)/T,t,-3,0); bn=int(2*sin(n*t*w)/T,t,0,1)+int(0*2*sin(n*t*w)/T,t,-3,0); ft=a0/2+symsum(an.*cos(n*t*w)+bn.*sin(n*t*w),n,1,10) t=-10:0.1:10;

ft=cos((pi*t)/2)/pi-cos((3*pi*t)/2)/(3*pi)+cos((5*pi*t)/2)/(5*pi)- cos((7*pi*t)/2)/(7*pi) +cos((9*pi*t)/2)/(9*pi)+sin(pi*t)/pi+sin((pi*t)/2)/pi+sin(3*pi*t)/(3*pi)+sin((3*pi*t)/2)/(3*pi)+sin(5*pi*t)/(5*pi)+sin((5*pi*t)/2)/(5*pi)+sin((7*pi*t)/2)/(7*pi)+sin((9*pi*t)/2)/(9*pi) + 1/4; plot(t,ft)

2.已知周期为T=4的三角波，在第一周期（-2

ft=0;T=4;w=2*pi/T;

a0=int(2*(t+1)/T,t,-2,0)+int(2*(1-t)/T,t,0,2);

an=int(2*(t+1)*cos(n*t*w)/T,t,-2,0)+int(2*(1-t)*cos(n*t*w)/T,t,0,2); bn=int(2*(t+1)*sin(n*t*w)/T,t,-2,0)+int(2*(1-t)*sin(n*t*w)/T,t,0,2); ft=a0/2+symsum(an.*cos(n*t*w)+bn.*sin(n*t*w),n,1,10)

32

t=-10:0.1:10;

ft=(8*cos((pi*t)/2))/pi^2+(8*cos((3*pi*t)/2))/(9*pi^2) + (8*cos((5*pi*t)/2))/(25*pi^2) + (8*cos((7*pi*t)/2))/(49*pi^2) + (8*cos((9*pi*t)/2))/(81*pi^2); plot(t,ft)

syms t n;

T1=4;w1=2*pi/T; a01=int(2/T1,t,0,1);

an1=int(2*cos(n*t*w1)/T1,t,0,1)+int(0*2*cos(n*t*w1)/T1,t,-3,0); bn1=int(2*sin(n*t*w1)/T1,t,0,1)+int(0*2*sin(n*t*w1)/T1,t,-3,0); Cn1=(an1-bn1*i)/2 for n=-20:20;

Cn1=sin((pi*n)/2)/(2*pi*n) - (sin((pi*n)/4)^2*i)/(pi*n); plot(w1*n,Cn1) hold on end

syms t n;

T2=4;w2=2*pi/T;

a02=int(2*(t+1)/T2,t,-2,0)+int(2*(1-t)/T2,t,0,2);

an2=int(2*(t+1)*cos(n*t*w2)/T2,t,-2,0)+int(2*(1-t)*cos(n*t*w2)/T2,t,0,2); bn2=int(2*(t+1)*sin(n*t*w2)/T2,t,-2,0)+int(2*(1-t)*sin(n*t*w2)/T2,t,0,2); Cn2=(an2-bn2*i)/2 for n=-20:20;

Cn2=(4*sin((pi*n)/2)^2)/(pi^2*n^2) - sin(pi*n)/(pi*n); plot(w2*n,Cn2) hold on end

33

本文来源：https://www.bwwdw.com/article/petw.html