信号与线性系统预备训练3——MATLAB软件在信号与系统中的应用初步

The Preparatory training3 of Signals and Linear Systems

对应教材：《信号与线性系统分析（第五版）》高等教育出版社，吴大正著

一、目的

1.熟悉和回顾MATLAB软件的基本操作和编程，为本课程翻转课堂和实验中基于MATLAB进行信号和系统的分析做好准备。
2.基于MATLAB，熟悉和回顾本课程所需的一些数学基础知识

二、任务安排及要求

1.根据所给出的各项内容要求，每位同学自主独立完成相关资料查阅、知识回顾和基于MATLAB的软件设计和调试，汇总分析相关结果并撰写报告（报告主要用于在课堂上做汇报展示，可以采用Word或者PPT形式，具体格式不限）。在完成任务的过程中，可以与其他同学商讨，或者向老师咨询，但应独立完成报告。
2.课堂上将随机抽取部分同学上台汇报展示（每位同学只展示其中部分内容），教师和其他同学针对该同学的汇报展示进行提问和讨论（教师提问有可能针对台下其他同学）。
3.熟练掌握MATLAB的基本操作和程序设计。

三、具体步骤要求

1、脚本与函数

要求：

（1）了解MATLAB的脚本和函数的概念，总结其联系和区别。举例说明对脚本文件和函数文件的命名有什么要求？
（2）任意编写一个MATLAB函数（函数功能自定），并利用脚本文件调用该函数，实现相应的函数功能。据此总结MATLAB函数的调用原理和方法。

解答：

（1）
联系：
脚本文件和函数文件都是m文件。
脚本是可以独立运行的，函数不能独立运行，需要被调用。
（2）
函数

%b、a分别为直接型的分子和分母系数
%C、B、A分别为并联型的常数项，分子项矩阵，分母项矩阵
function [C,B,A]=tf2par(b,a)
[r,p,C]=residuez(b,a);
N=length(r);            
A=zeros(ceil(N/2),3);
B=A;                    
k=1;                   
if mod(N,2)==1          
for i=1:2:N-1
[B(k,:),A(k,:)]=residuez([r(i),r(i+1)],[p(i),p(i+1)],0);
k=k+1;
end
B(end,:)=[r(end),0,0];
A(end,:)=[1,-p(end),0];
else                    
for i=1:2:N
[B(k,:),A(k,:)]=residuez([r(i),r(i+1)],[p(i),p(i+1)],0);
k=k+1;
end
end
end

调用

syms z;
A=[1 2 -3 4 5 -6 7];
B=[7 -6 5 4 -3 2 1];
[S,G]=tf2sos(B,A)
[c,b,a]=tf2par(B,A)

调用格式：
[输出参数1，输出参数2，…]=函数名（输入参数1，输入参数2，…）
函数调用可以嵌套，一个函数可以调用别的函数，甚至调用它自己。

2、符号运算和数值运算

要求：

符号运算和数值运算是MATLAB中两种不同的数学运算模式，其用途也不同。符号运算一般用于理论推导，得到的是运算结果的解析表达式；数字运算一般用于数值求解和图形化展示，得到的运算结果是数据（数组、矩阵），可以直接绘制成图表的形式。
（1）进一步了解符号运算与数值运算的区别，掌握两种不同运算模式的基本编程方法。
（2）了解符号运算函数dsolve（求解微分方程）、fourier（傅里叶变换）、ifourier（傅里叶反变换）、laplace（拉普拉斯变换）、ilaplace（拉普拉斯反变换）、ztrans（Z变换）、iztrans（Z反变换）等的功能，分别举例说明这些函数的用法。
（3）根据傅里叶变换的定义，采用数值运算的方法计算某个函数（具体函数自定）的傅里叶变换。若自己独立编程完成有困难，可以咨询老师或者到网络上搜索相关示例，但应该理解其运算原理和程序设计原理。

解答：

（1）
创建符号变量的基本方法——利用syms命令
例如：

syms x a b
f=x^2+a*x+b;
fx=diff(f,x)

运行结果：

fx = a + 2*x

数值运算为直接赋值
例如：

a=1;
b=2;
c=a+b

运行结果：

c = 3

（2）
dsolve：用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。
句法：

S = dsolve(eqn)
S = dsolve(eqn,cond)
[y1,...,yN] = dsolve(___)

句法解释：

S＝dsolve（eqn）求解微分方程eqn，其中eqn是符号方程。使用diff和==表示微分方程。例如，diff（y，x）==y表示等式dy/dx=y。通过将方程指定为微分方程的向量来求解微分方程组。
例子：
解方程：
$\frac{\mathrm{d}y}{\mathrm{d}t}=ay$

syms a y(t)
eqn = diff(y,t) == a*y;
dsolve(eqn)

ans =
C1*exp(a*t)

解高阶方程：
$\frac{\mathrm{d^{2}}y}{\mathrm{d}t^{2}}=ay$

syms y(t) a
eqn = diff(y,t,2) == a*y;
ySol(t) = dsolve(eqn)

ySol(t) =
C2*exp(-a^(1/2)*t) + C3*exp(a^(1/2)*t)

S＝dsolve（eqn，cond）用初始或边界条件cond求解方程。
例子：
解如下方程：
$\frac{\mathrm{d}y}{\mathrm{d}t}=ay$
且已知条件：
$y （ 0 ） = 5$

syms y(t) a
eqn = diff(y,t) == a*y;
cond = y(0) == 5;
ySol(t) = dsolve(eqn,cond)

ySol(t) =
5*exp(a*t)

解高阶方程：
$\frac{\mathrm{d^{2}}y}{\mathrm{d}t^{2}}=a^{2}y$
且条件：
$y (0) = b, y^{'} (0) = 1$

syms y(t) a b
eqn = diff(y,t,2) == a^2*y;
Dy = diff(y,t);
cond = [y(0)==b, Dy(0)==1];
ySol(t) = dsolve(eqn,cond)

ySol(t) =
(exp(a*t)*(a*b + 1))/(2*a) + (exp(-a*t)*(a*b - 1))/(2*a)

[y1，…，yN]=dsolve（__）将解分配给变量y1，…，yN。

例子：
解方程组：
$\begin{cases} \begin{aligned} \frac{\mathrm{d}y}{\mathrm{d}t}=z\\ \frac{\mathrm{d}z}{\mathrm{d}t}=-y \end{aligned} \end{cases}$

syms y(t) z(t)
eqns = [diff(y,t)==z, diff(z,t)==-y];
[ySol(t) zSol(t)] = dsolve(eqns)

ySol(t) =
C2*cos(t) + C1*sin(t)
zSol(t) =
C1*cos(t) - C2*sin(t)

fourier：傅里叶变换——将满足一定条件的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合。
F=fourier（f）是符号表达式的傅立叶变换或具有默认自变量x的函数f。默认情况下，结果F是w的函数。如果f=f（w），则F作为变量v的函数返回F＝F（v）。
F=fourier（f，v）返回F作为变量v的函数而不是默认变量w。
F=fourier（f，u，v）将F视为变量u的函数而不是默认变量x。

syms t v w x f(x)fourier(1/t)   returns   -pi*sign(w)*1ifourier(exp(-x^2),x,t)   returns   pi^(1/2)*exp(-t^2/4)fourier(exp(-t)*heaviside(t),v)   returns   1/(1+v*1i)fourier(diff(f(x)),x,w)   returns   w*fourier(f(x),x,w)*1i

ifourier：傅里叶逆变换
f＝ifourier（F）是符号的傅立叶逆变换。带有默认自变量w的表达式或函数F。结果f默认是x的函数。如果F＝F（x），则f作为变量t的函数返回，f＝f（t）。
f=ifourier（F，u）返回f作为变量u的函数,而不是默认变量x。
f=ifourier（F，v，u）将F视为变量v的函数而不是默认变量w。
例如：

syms t u v w f(x)ifourier(w*exp(-3*w)*heaviside(w))  returns  1/(2*pi*(-3+x*1i)^2)ifourier(1/(1 + w^2),u)   returns   exp(-abs(u))/2ifourier(v/(1 + w^2),v,u)   returns   -(dirac(1,u)*1i)/(w^2+1)ifourier(fourier(f(x),x,w),w,x)   returns   f(x)

laplace：以实现信号f(t)的拉普拉斯变换
L=laplace (f) 是符号函数f的拉普拉斯变换，默认返回函数L是关于s的函数。
L=laplace (f,t) 是符号函数f的拉普拉斯变换，返回函数L是关于t的函数。
L=laplace (f,w,z) 是关于w的函数f的拉普拉斯变换，返回函数L是关于z的函数。
函数调用实例

syms a s t w xlaplace(t^5)                 returns   120/s^6laplace(exp(a*s))             returns   1/(t-a)laplace(sin(w*x),t)            returns   w/(t^2+w^2)laplace(cos(x*w),w,t)          returns   t/(t^2+x^2)laplace(x^sym(3/2),t)          returns   3/4*pi^(1/2)/t^(5/2)laplace(diff(sym('F(t)')))        returns   laplace(F(t),t,s)*s-F(0)

ilaplace：可以实现信号F(s)的拉普拉斯逆变换
F = ilaplace(L) 是函数L(s)的拉普拉斯逆变换，默认返回函数F是关于t的函数。
F = ilaplace(L,y) 是函数L(s)的拉普拉斯逆变换，返回函数F是关于y的函数。
F = ilaplace(L,y,x) 返回函数F是关于x的函数。
函数调用实例

 syms s t w x yilaplace(1/(s-1))                  returns   exp(t)ilaplace(1/(t^2+1))                returns   sin(x)ilaplace(t^(-sym(5/2)),x)            returns   4/3/pi^(1/2)*x^(3/2)ilaplace(y/(y^2 + w^2),y,x)          returns   cos(w*x)ilaplace(sym('laplace(F(x),x,s)'),s,x)   returns   F(x)

ztrans:实现信号f(k)的z变换
F = ztrans (f) 是函数f(n)的z变换，默认返回函数F是关于z的函数。
F = ztrans (f,w) 是函数f(n)的z变换，返回函数F是关于w的函数。
F = ztrans (f,k,w) 是函数f(k)的z变换，返回函数F是关于w的函数。
函数调用范例

syms k n w z
ztrans(2^n)                returns  z/(z-2)
ztrans(sin(k*n),w)          returns  sin(k)*w/(1-2*w*cos(k)+w^2)
ztrans(cos(n*k),k,z)         returns  z*(-cos(n)+z)/(-2*z*cos(n)+z^2+1)
ztrans(cos(n*k),n,w)         returns  w*(-cos(k)+w)/(-2*w*cos(k)+w^2+1)
ztrans(sym('f(n+1)'))         returns  z*ztrans(f(n),n,z)-f(0)*z

iztrans:实现信号F(z)的逆z变换
f = iztrans (F) 是函数F(z)的逆z变换，默认返回函数f是关于n的函数。
f = iztrans (F,k) 是函数F(z)的逆z变换，返回函数f是关于k的函数。
f = iztrans (F,w,k) 是函数F(w)的逆z变换，返回函数f是关于k的函数。
函数调用范例：

syms z x k f(n)
iztrans(z/(z-2))        returns   2^n
iztrans(sin(1/n))  returns  -(1i^(k-1)*((-1)^k-1))/(2*factorial(k))
iztrans(exp(x/z),z,k)    returns   x^k/k
iztrans(ztrans(f(n),n,z),z,k)  returns  f(k)

（3）
对函数x=sin(20πt)+sin(80πt)做傅里叶变换

t=0:1/100:10-1/100;
x=sin(2*pi*10*t)+sin(2*pi*40*t);
y=fft(x);
m=abs(y);
f=(0:numel(y)-1)'*100/numel(y);
plot(f,m);
grid on;
ylabel('幅度');
xlabel('频率');

在这里插入图片描述

3、运算结果的可视化函数调用范例

要求：

将运算结果表达成更为直观的图形化形式，是MATLAB的一个重要功能，便于直观地理解一些比较抽象的概念和理论分析结果。一般来说，图形化表示应基于数值运算的结果，但对于符号运算结果，MATLAB也定义了一些直接基于解析表达式来绘图的函数。
运算结果图形化表示时需要特别注意的一个问题，是绘图坐标轴范围的合理选取问题：图形化表示的坐标范围总是有限的，应该在有限的绘图区间内，充分展现所关注的运算结果的特征。绘图坐标范围太小或者区间选择不合理，不能完整展现运算结果的特征；坐标范围太大，则无关图形太多，特征图形太小，其特性展现不突出。
（1）基于数值运算结果的常用可视化函数。了解plot（二维曲线）、plot3（三维曲线）、surf（三维阴影曲面）、mesh（三维网格曲面）、stem（二维离散序列）、stem3（三维离散序列）等常用图形化函数的功能和用法，并举例展示。
（2）基于符号运算结果的常用可视化函数。了解ezplot（二维曲线）、ezplot3（三维曲线）、ezsurf（三维阴影曲面）、ezmesh（三维网格曲面）等基于解析表达式的图形化函数的功能和用法，并利用这些函数分别绘制dsolve（求解微分方程）、fourier（傅里叶变换）、laplace（拉普拉斯变换）等符号函数运算结果的图形，注意根据需要选择适当的函数、维数和绘图坐标轴范围。
（3）给定一个函数（自定），分别用符号运算和数字运算的方式求其傅里叶变换，然后分别基于数值运算和符号运算结果绘图，注意选择适当的数值计算范围和绘图坐标轴范围。

解答：

（1）
plot

x = 0:pi/10:2*pi;
y1 = sin(x);
y2 = sin(x-0.25);
y3 = sin(x-0.5);figure
plot(x,y1,'g',x,y2,'b--o',x,y3,'c*')

在这里插入图片描述
plot3

t = 0:pi/50:10*pi;
st = sin(t);
ct = cos(t);figure
plot3(st,ct,t)

在这里插入图片描述
surf

[X,Y] = meshgrid(1:0.5:10,1:20);
Z = sin(X) + cos(Y);
surf(X,Y,Z)

在这里插入图片描述
mesh

[X,Y] = meshgrid(-8:.5:8);
R = sqrt(X.^2 + Y.^2) + eps;
Z = sin(R)./R;
C = gradient(Z);figure
mesh(X,Y,Z,C)

在这里插入图片描述
stem

figure
X = linspace(0,2*pi,25)';
Y = (cos(2*X));
stem(X,Y,'LineStyle','-.',...'MarkerFaceColor','red',...'MarkerEdgeColor','green')

在这里插入图片描述
stem3

figure
theta = linspace(0,2*pi);
X = cos(theta);
Y = sin(theta);
Z = theta;
stem3(X,Y,Z,':*m')

在这里插入图片描述

（2）
ezplot

ezplot('x^2-y^4')

在这里插入图片描述
ezplot3

ezplot3('sin(t)','cos(t)','t',[0,6*pi])

在这里插入图片描述

ezsurf

ezsurf('sqrt(1-x^2-y^2)')

在这里插入图片描述
ezmesh

fh = @(x,y) x.*exp(-x.^2-y.^2);
ezmesh(fh,40)

在这里插入图片描述

（3）

暂无，欢迎贡献分享

4、音频信号采集和保存、回放

要求：

MATLAB可以调用电脑的声卡进行音频信号的采集和回放。自行查阅相关资料，了解基于audiorecorder对象的音频采集、保存（含保存为音频文件）、音频文件读取、回放等功能的程序实现方法，并举例说明。

解答：

music=audiorecorder(8000,16,2);
%创建一个保存音频信息的对象，它包含采样率，时间和录制的音频信息等等。
%44100表示采样为44100Hz（可改为8000, 11025, 22050等，
%此数值越大，录入的声音质量越好，相应需要的存储空间越大）
%16为用16bits存储，2为两通道即立体声（也可以改为1即单声道）。recordblocking(music,5);
%开始录制，此时对着麦克风说话即可,录制时间为5秒。
play(music);MyRecording=getaudiodata(music);
%得到以n*2列数字矩阵存储的刚录制的音频信号。
plot(MyRecording);filename='myspeech.wav';
audiowrite(filename,MyRecording,8000);
%MyRecording表示要存入的波形矩阵，
%8000表采样率，'myspeech'为存储的文件名
[y, fs]=audioread('myspeech.wav');

在这里插入图片描述

5、帮助和错误提示

要求：

（1）MATLAB内部定义了大量的功能函数供用户使用，当用户不清楚函数的功能和用法时，可以在命令行窗口用“help 函数名”查看该函数的帮助信息，并通过其中的链接了解该函数的更多细节（英文）。请举例说明help命令的使用。
（2）当执行的MATLAB命令（直接在命令行运行）或者脚本存在错误时，MATLAB会给出相应的错误提示，以便用户排除错误。请查阅网络资料，了解MATLAB中常见的错误提示及其产生的原因，并通过例子说明其中的一些错误。