文章目录
一、实验目的：
二、实验内容及要求：
1. 离散时间信号的频域分析：
2. 线性时不变离散时间系统的频域分析：
三、实验结果及问题回答：
1. 离散时间信号的频域分析：
2. 线性时不变离散时间系统的频域分析：

一、实验目的：
熟悉Mat lab中的信号处理工具箱及相关命令，掌握用Mat lab验证离散时间傅里叶变换相关性质的方法；掌握用Mat lab分析离散时间系统频域特性的方法并能绘制幅频特性和相频特性图。
二、实验内容及要求：
1. 离散时间信号的频域分析：
（1）修改程序P3.1，在范围 内计算如下序列的离散时间傅里叶变换并绘制幅频特性和相频特性图。



（2）修改程序P3.2，D=8，验证离散时间傅里叶变换的时移特性并绘图。

（3）修改程序P3.3，w = 0.2π ，验证离散时间傅里叶变换的频移特性并绘图。
2. 线性时不变离散时间系统的频域分析：
（1）修改程序P3.1，计算并画出如下传输函数的频率响应，它是哪种类型的滤波器？



（2）修改程序P4.2，计算并画出一个长度为6的滑动平均滤波器的增益响应。

（3）用Matlab画出如下滤波器的零极点图并分析稳定性。


三、实验结果及问题回答：
1. 离散时间信号的频域分析：
（1）实验结果：
% Program P3_1
% Evaluation of the DTFT 
clf;
% Compute the frequency samples of the DTFT
w = 0*pi:8*pi/511:pi;
num = [0.7 -0.5 0.3 1];den = [1 0.3 -0.5 0.7];
h = freqz(num, den, w);
% Plot the DTFT
subplot(2,1,1)
plot(w/pi,real(h));grid
title('Real part of H(e^{j\omega})')
xlabel('\omega /\pi');
ylabel('Amplitude');
subplot(2,1,2)
plot(w/pi,imag(h));grid
title('Imaginary part of H(e^{j\omega})')
xlabel('\omega /\pi');
ylabel('Amplitude');
pause
subplot(2,1,1)
plot(w/pi,abs(h));grid
title('Magnitude Spectrum |H(e^{j\omega})|')
xlabel('\omega /\pi');
ylabel('Amplitude');
subplot(2,1,2)
plot(w/pi,angle(h));grid
title('Phase Spectrum arg[H(e^{j\omega})]')
xlabel('\omega /\pi');
ylabel('Phase in radians');



（2）实验结果：
% Program P3_2
% Time-Shifting Properties of DTFT
clf;
w = -pi:2*pi/255:pi; wo = 0.4*pi; D = 8;
num = [1 2 3 4 5 6 7 8 9];
h1 = freqz(num, 1, w);
h2 = freqz([zeros(1,D) num], 1, w);
subplot(2,2,1)
plot(w/pi,abs(h1));grid
title('Magnitude Spectrum of Original Sequence')
subplot(2,2,2)
plot(w/pi,abs(h2));grid
title('Magnitude Spectrum of Time-Shifted Sequence')
subplot(2,2,3)
plot(w/pi,angle(h1));grid
title('Phase Spectrum of Original Sequence')
subplot(2,2,4)
plot(w/pi,angle(h2));grid
title('Phase Spectrum of Time-Shifted Sequence')



（3）实验结果：
% Program P3_3
% Frequency-Shifting Properties of DTFT
clf;
w = -pi:2*pi/255:pi; wo = 0.2*pi;
num1 = [1 3 5 7 9 11 13 15 17];
L = length(num1);
h1 = freqz(num1, 1, w);
n = 0:L-1;
num2 = exp(wo*i*n).*num1;
h2 = freqz(num2, 1, w);
subplot(2,2,1)
plot(w/pi,abs(h1));grid
title('Magnitude Spectrum of Original Sequence')
subplot(2,2,2)
plot(w/pi,abs(h2));grid
title('Magnitude Spectrum of Frequency-Shifted Sequence')
subplot(2,2,3)
plot(w/pi,angle(h1));grid
title('Phase Spectrum of Original Sequence')
subplot(2,2,4)
plot(w/pi,angle(h2));grid
title('Phase Spectrum of Frequency-Shifted Sequence')


2. 线性时不变离散时间系统的频域分析：
（1）实验结果：

% Program P3_1

% Evaluation of the DTFT

clf;

% Compute the frequency samples of the DTFT

w = 0pi:8pi/511:1*pi;

num = [0.15 0 -0.15];den = [1 -0.5 0.7];

h = freqz(num, den, w);

% Plot the DTFT

subplot(2,1,1)

plot(w/pi,real(h));grid

title(‘Real part of H(e^{j\omega})’)

xlabel(’\omega /\pi’);

ylabel(‘Amplitude’);

subplot(2,1,2)

plot(w/pi,imag(h));grid

title(‘Imaginary part of H(e^{j\omega})’)

xlabel(’\omega /\pi’);

ylabel(‘Amplitude’);

pause

subplot(2,1,1)

plot(w/pi,abs(h));grid

title(‘Magnitude Spectrum |H(e^{j\omega})|’)

xlabel(’\omega /\pi’);

ylabel(‘Amplitude’);

subplot(2,1,2)

plot(w/pi,angle(h));grid

title(‘Phase Spectrum arg[H(e^{j\omega})]’)

xlabel(’\omega /\pi’);

ylabel(‘Phase in radians’);


（2）实验结果：
% Program P4_2
% Gain Response of a Moving Average Lowpass Filter
clf;
M = 6;
num = ones(1,M)/M;
[g,w] = gain(num,1);
plot(w/pi,g);grid
axis([0 1 -50 0.5])
xlabel('\omega /\pi');ylabel('Gain in dB');
title(['M = ', num2str(M)])

function [g,w] = gain(num,den)
% Computes the gain function in dB of a 
% transfer function at 256 equally spaced points 
% on the top half of the unit circle
% Numerator coefficients are in vector num
% Denominator coefficients are in vector den
% Frequency values are returned in vector w
% Gain values are returned in vector g
w = 0:pi/255:pi;
h = freqz(num,den,w);
g = 20*log10(abs(h));



（3）实验结果：
clc;

clear;

num = [ 1 0 0];

den = [1 -1.848 0.85];

zplane(num,den)

信号与系统的分析方法有两种：时域分析方法和频域分析方法。
在连续时间信号与系统中，信号一般用连续变量时间t的函数表示，系统用微分方程描述，其频域分析方法是拉普拉斯变换和傅立叶变换。在时域离散信号与系统中，信号用序列表示，其自变量仅取整数，非整数时无定义，系统则用差分方程描述，频域分析方法是Z变换和序列傅立叶变换法。
Z变换在离散时间系统中的作用就如同拉普拉斯变换在连续时间系统中的作用一样，它把描述离散系统的差分方程转化为简单的代数方程，使其求解大大简化。因此，对求解离散时间系统而言，Z变换是一个极重要的数学工具。
2．2序列的傅立叶变换(离散时间傅立叶变换)
一、序列傅立叶变换：
正变换：DTFT[x(n)]=(2.2.1)
反变换：DTFT-1
式(2.2.1)级数收敛条件为
||=   (2.2.2)
上式称为x(n)绝对可和。这也是DTFT存在的充分必要条件。当遇到一些绝对不可和的序列，例如周期序列，其DTFT可用冲激函数的形式表示出来。
二、序列傅立叶变换的基本性质：
1、DTFT的周期性
，是频率w的周期函数，周期为2p。
∵ = 。
问题1：设x(n)=RN(n)，求x(n)的DTFT。
==
==
设N为4，画出幅度与相位曲线。
2、线性
设=DTFT[x1(n)]，=DTFT[x2(n)]，
则：DTFT[a x1(n)+b x2(n)]
=  = a+b
3、序列的移位和频移
设 = DTFT[x(n)]，
则：DTFT[x(n-n0)] = 
= 
DTFT[x(n)]  = 
=  = 
4、DTFT的对称性
共轭对称序列的定义：设序列满足下式
则称为共轭对称序列。
共轭对称序列的性质：
共轭对称序列的实部是偶函数，虚部是奇函数
证明：=+j(实部加虚部)
∵
∴+j=-j
∴=(偶函数)
∴=－(奇函数)
一般情况下，共轭对称序列用表示：
共轭反对称序列的定义：设序列满足下式
则称为共轭反对称序列。
共轭反对称序列的性质：
共轭反对称序列的实部是奇函数，虚部是偶函数
证明：=+j(实部加虚部)
∵
∴+j=+j
∴=(奇函数)
∴=(偶函数)
一般情况下，用来表示
一个序列可用共轭对称序列与共轭反对称序列之和表示。即：
x(n)= +     (2.2.16)
问题1： =？
=+
=－
∴= －    (2.2.17)
=(+)
=(－)
对于频域函数，也可分解成共轭对称分量和共轭反对称分量之和：
式中，是共轭对称分量，是共轭反对称分量，它们满足：
=，=
且：
：共轭对称分量，它的实部是偶函数，虚部是奇函数；：共轭反对称分量，它的实部是奇函数，虚部是偶函数。
下面研究DTFT的对称性，按下面两部分进行分析
a)将序列x(n)分成实部与虚部，即：
=+j(、都是实数序列)
则：
式中：=DTFT[]=，
=DTFT[j]=j。
结论：序列分成实部与虚部两部分，实部对应于中的，虚部和j一起对应于中的。
b)将序列分成共轭对称部分和共轭反对称部分，x(n)= + 
∵=(+)
=(－)
将上面两式分别进行DTFT，得到：
DTFT[]=(+)=Re[]=
DTFT[]=()=jIm[]=j
∴=+j
x(n)= + 
结论：序列的共轭对称部分对应于的实部，而序列的共轭反对称部分对应于的虚部加j。
应用：利用DTFT的对称性讨论当h(n)是实序列时，其DTFT的特性。
∵h(n)是实序列，所以它所对应的DTFT：=，具有共轭对称性，的实部偶对称，虚部奇对称。
5、时域卷积定理
设y(n)=x(n)*h(n)
则：=×=  (2.2.32)
证明：y(n)= x(n)*h(n)=
=DTFT[y(n)]
= =
=
=
=
6、频域卷积定理
设y(n) = x(n) h(n)
则=*=
=
证明：==
=
=
=
=*
7、Parseval(帕斯维尔)(帕塞瓦尔)定理
= (2.2.34)
证明：
==
= 
=
=
2．5  Z变换的定义与收敛域
一、Z变换的定义
若序列为x(n)，则幂级数
(2.5.1)
称为序列x(n)的Z变换，也称为双边Z变换。式中z为复变量，它所在的复平面称为z平面。亦可将x(n)的Z变换表示为
ZT[x(n)] = X(z)
二、Z变换的收敛域
我们知道，是一幂级数，只有收敛时Z变换才有意义。X(z)收敛的条件是：
      (2.5.3)
X(z)能够收敛的z取值集合称为X(z)的收敛域。
一般收敛域用环状域表示。即：
∴Z变换的公式
   (2.5.1)
常见的Z变换是一个有理函数，表示为：
分子多项式的根是的零点，分母多项式的根是的极点。在极点处Z变换不存在。因此收敛域中没有极点，收敛域总是用极点限定其边界。
1、有限长序列Z变换的收敛域
有限长序列是指在有限区间n1≤n≤n2之间序列具有非零的有限值，在此区间外，序列值皆为零。有限长序列Z变换为，所以收敛域为
0
如n1≥0，收敛域为0
如n2≤0，收敛域为0≤|z|
2、右边序列Z变换的收敛域
右边序列是指在n≥n1时，x(n)有值，在n＜n1时，x(n)=0。其Z变换为
此式右端第一项为有限长序列的Z变换，它的收敛域为0≤|z|。综合此两项，只有两项都收敛时级数才收敛。所以右边序列Z变换的收敛域为。
因果序列是最重要的一种右边序列，即n1=0的右边序列。收敛域为(也可以写成)，所以，|z|=∞处Z变换收敛是因果序列的特征。
3、左边序列Z变换的收敛域
左边序列是指在n≤n2时，x(n)有值，n＞n2时，x(n)=0。其Z变换为
此式第二项是有限长序列的Z变换，收敛域为0
4、双边序列Z变换的收敛域
这类序列是指n为任意值时x(n)皆有值的序列。双边序列的收敛域为
问题1：求序列x(n)= RN(n)的Z变换及收敛域，并画出收敛域。
解：X(z)==。因为这是有限长序列，所以收敛域为0
思考：RN(n)的DTFT存在吗？
问题2：x(n)=anu(n)，求其Z变换及收敛域，并画出收敛域。
解：这是右边序列，且是因果序列，其Z变换为X(z)=。收敛域为(或写成)
思考：anu(n)的DTFT存在吗？
问题3：x(n)=-anu(-n-1)，求其Z变换及收敛域，并画出收敛域。
解：这是一个左边序列。其Z变换为
，
收敛域为0≤|z|
思考：-anu(-n-1)的DTFT存在吗？
结论：当Z变换的收敛域中包含单位圆时，用Z变换可求出DTFT。
=(2.5.4)
上式称为单位圆上的Z变换就是离散时间傅立叶变换。
回顾：观察零极点。
结论：零点可以在复平面的任意处，但极点在收敛域的边缘或收敛域的外面。
2.5.3  Z反变换
已知序列的Z变换及其收敛域，求序列称为Z反变换。表示为x(n)=ZT-1[X(z)]
         
其中，c是X(z)收敛域中一条逆时针的闭合曲线。
求Z反变换的方法通常有三种：围线积分法(留数法)、部分分式展开法和长除法。
一、围线积分法(留数法)
直接计算围线积分比较麻烦，一般都采用留数定理来求解。按留数定理，若函数F(z)=X(z)zn-1在围线c上连续，在c以内有K个极点zk，则有
(2.5.6)
设zr是X(z)zn-1的单极点，则根据留数定理：
如果zk是L阶极点，则根据留数定理，
(2.5.8)
(2.5.8)表明，对于L阶极点，需要求L-1次导数，这是比较麻烦的。如果c内有多阶极点，而c外没有多阶极点时，可根据留数辅助定理改求c外所有极点之和，使问题简单化。
若函数F(z)=X(z)zn-1在围线c上连续，在c以内有K个极点zk，而在c以外有M个极点zm，(K，M为有限值)。现在c内有多阶极点，而c外没有多阶极点，根据留数辅助定理改求c外所有极点之和。得：
(2.5.9)(2.5.9)应用条件是X(z)zn-1在z=∞有两阶或二阶以上零点，即要分母多项式z的阶次比分子多项式z的阶次高二阶或二阶以上。
问题1：已知X(z)=z2/[(4-z)(z-1/4)]，1/4
求Z反变换。
解：c，c为X(z)的收敛域
内的闭合围线，画出收敛域及c。
X(z)zn-1=。现在来看极点在围线c内部及外部的分布情况及极点阶数。
当时，函数在围线c内只有z=1/4处一个一阶极点，
=，
当时，函数在围线外部只有一个一阶极点z=4，而在围线的内部则有z=1/4处一阶极点及z=0处一(n+1)阶极点，所以采用围线外部的极点较方便。
=，
∴
问题2：已知X(z)=z2/[(4-z)(z-1/4)]，|z|>4，
求Z反变换。
解：c，c为X(z)的收敛域
内的闭合围线。
X(z)zn-1=。现在来看在围c内部及外部的分布情况及极点阶数。
当时，
函数在围线c内z=1/4处有一个一阶极点，z=4处有一个一阶极点，
+
=，
当n=-1时，x(n)=0，∴x(n)= ，
当时，函数在围线外部没有一个极点，所以采用围线外部的极点较方便。由于围线外部没有一个极点，∴x(n)=0。
∴x(n)=()u(n)
二、部分分式展开法
对于大多数单极点的序列，常常用这种部分分式展开法求Z反变换。
X(z)=B(z)/A(z)= X1(z)+ X2(z)+…+ XK(z)，则
=ZT-1[X1(z)]+ZT-1[ X2(z)]+…+ZT-1[XK(z)]
ZT-1[X1(z)]、ZT-1[ X2(z)]、…ZT-1[XK(z)]可从Z变换表中直接查表得出
问题1：设X(z)=z2/[(z-2)(z-0.5)]，|z|>2，
求Z反变换。
解：X(z) =z2/[(z-2)(z-0.5)]
A1=，A2=
∴，
∵收敛域为|z|>2，∴x(n)=
三、幂级数展开法
因为的Z变换定义为z-1的幂级数，即
所以只要在给定得收敛域内，把X(z)展成幂级数，则级数的系数就是序列。
当X(z)的收敛域为|z|>Rx-时，则必为因果序列，此时应将X(z)展成z的负幂级数，为此，X(z)的分子分母应按z的降幂排列；
当X(z)的收敛域为|z|必为左边序列，X(z)的分子分母应按z的升幂排列；
问题1：已知，|z|>3
解：因为收敛域|z|>3，所以这是因果序列，因此，X(z)分子分母按z的降幂排列。
进行长除
2.5.4  Z变换的基本性质和定理
一、线性
线性就是要满足比例性和可加性。若
X(z) =ZT [x(n) ]，
Y(z) = ZT [y(n) ]，
则ZT[ax(n)+by(n)]=aX(z)+b Y(z)，
，。
二、序列的移位
若X(z) = ZT [x(n) ]，
则有ZT [x(n-m) ] =z-mX(z)，
三、乘以指数序列
若X(z) =ZT [x(n) ]，
则ZT [anx(n) ]=X()，
四、序列乘以n
若X(z) =ZT [x(n) ]，
则ZT [n x(n) ]=－z，
五、复序列取共扼
一个复序列x(n)的共扼序列为x*(n)
若ZT [x(n) ] =X(z)，
则ZT [x*(n) ] =X*(z*)，
六、翻转序列
若ZT [x(n) ] =X(z)，
则ZT [x(－n) ] =X()，
七、(因果序列)初值定理
对于因果序列x(n)，即x(n)=0，n<0，ZT[x(n) ] =X(z)有
八、(因果序列)终值定理
设x(n)为因果序列，且X(z) = ZT [x(n) ]的极点处于单位圆|z|=1以内(单位圆上最多在z=1处可有一阶极点)，则
九、序列的卷积和(时域卷积和定理)
设y(n)为x(n)与h(n)的卷积和
y(n)= x(n)*h(n)=
X(z) =ZT [x(n) ]，
H(z) =ZT [h(n) ]，
则Y(z) =ZT [y(n) ]= X(z) H(z)，
十、序列相乘(z域卷积定理)
若y(n)= x(n)·h(n)，且
X(z) =ZT [x(n) ]，
H(z) =ZT [h(n) ]，
则Y(z) =ZT [y(n) ]=ZT[x(n)·h(n)]
=，
其中c是v平面上，与H(v)的公共收敛域内环绕原点的一条反时针旋转的单封闭围线。
v平面收敛域为
或Y(z) =ZT [y(n) ]=ZT[x(n)·h(n)]
=，
其中c是v平面上，与X(v)的公共收敛域内环绕原点的一条反时针旋转的单封闭围线。
v平面收敛域为
十一、帕斯维尔(Parseval)定理
若X(z) =ZT [x(n) ]，
H(z) =ZT [h(n) ]，
且
则
v平面上，c所在的收敛域为
证明：Y(z) =ZT[x(n)·h*(n)]
=
=，
因为，所以z=1在收敛域中。令z=1代入上式，
=
v平面上，c所在的收敛域为
如果X(z)，H(z)在单位圆上都收敛，则c可取为单位圆，即，则
如果h(n)=x(n)，则进一步有。
2.5.5利用Z变换解差分方程
在第一章中介绍了差分方程的递推解法，下面介绍Z变换解法。这种方法将差分方程变成了代数方程，使求解过程简单
设N阶线性常系数差分方程为
     (2.5.30)
一、求及
对(2.5.30)求双边Z变换：
=
=/=，h(n)=ZT-1[]
=，y(n)=ZT-1[]
2.6离散系统的系统函数，系统的频率响应
信号和系统的频率特性一般用序列的傅立叶变换和Z变换进行分析。
一、传输函数与系统函数
设系统初始状态为零，输出端对输入为单位抽样序列d(n)的响应，称为系统的单位抽样响应h(n)。对h(n)进行傅立叶变换得到：=，一般称为为系统的传输函数，它表征系统的频率特性。
将h(n)进行Z变换，得到，一般称H(z)为系统的系统函数，它表征了系统的复频域特性。
如已知系统的N阶线性常系数差分方程，进行双边Z变换，得到系统函数的一般表示式：
如果的收敛域包含单位圆|z|=1则，与的关系：=。即单位圆上的系统函数就是系统的传输函数
二、用系统函数的极点分布分析系统的因果性和稳定性
因果(可实现)系统其单位脉冲响应h(n)一定满足：当n<0时，h(n)=0，那么其系统函数的收敛域一定包含¥点。
系统稳定要求，对照ZT定义，系统稳定要求收敛域包含单位圆。
所以系统因果且稳定，收敛域包含¥点和单位圆，那么收敛域表示为：r
问题1：一个因果系统的系统函数为=，其中a为实数，问：a在哪些范围内才能使系统稳定？
解：因为系统因果，所以收敛域为|a|
三、利用系统的零极点分布分析系统的频率特性
=
将上式因式分解，得到：
=A
式中，是的零点，是其极点。A参数影响传输函数的幅度大小，影响系统特性的是零点和极点的分布。下面采用几何方法研究系统零极点分布对系统频率特性的影响。
=A
设系统稳定，将z=代入，得：
=A
在z平面上，-用一根由零点指向单位圆上点B的向量表示。同样-用由极点指向点B的向量表示，如图2.6.2。
将向量用极坐标表示：=，=，得到：
=A=
= |A|   (2.6.8)
=(N=M)(2.6.9)
当频率w从零变化到2p时，这些向量的终点B沿单位圆逆时针旋转一周，按照(2.6.8)和(2.6.9)分别估算出系统的幅度特性和相位特性。
按照(2.6.8)知道零极点的分布后，可以很容易地确定零极点位置对系统特性的影响。当B点转到极点附近时，极点矢量长度最短因而幅度特性可能出现峰值，且极点愈靠近单位圆，极点矢量长度愈短，峰值愈高愈尖锐。如果极点在单位圆上，则幅度特性为¥，系统不稳定。对于零点，情况相反，当B点转到零点附近，零点矢量长度变短，幅度特性将出现谷值，零点愈靠近单位圆，谷值愈接近零。当零点处在单位圆上时，谷值为零。总结以上结论：极点位置主要影响频响的峰值位置及尖锐程度，零点位置主要影响频响的谷点位置及形状。1．5模拟信号数字处理方法
数字信号处理技术优于模拟信号处理技术，故人们将模拟信号数字化，即经过采样、量化编码最终形成数字信号。
连续时间信号变为离散时间信号是由“采样”这一过程完成的。采样是将模拟信号数字化的第一个环节。它是利用周期性抽样脉冲序列(常用p(t)表示)从连续信号中抽取一系列的离散值来得到抽样信号的。如下图，根据每个脉冲宽度的不同，可将抽样分为两种:
理想抽样实际抽样
我们要研究的是，信号被抽样后其频谱将会有什么变化？在什么条件下，可从抽样信号中不失真地恢复原来信号xa(t)?
设：xa(t)的傅立叶变换为：，抽样脉冲序列p(t)的傅立叶变换为：，抽样信号的傅立叶变换为：，
∵= xa(t)p(t)，∴。
由上图得，抽样脉冲序列p(t)的周期为T，则抽样频率。则周期信号p(t)的傅立叶变换，其中。
的傅立叶变换为：
一、理想抽样
p(t)=
=
=
∴
=
= xa(t) p(t)，p(t)=
= xa(t) p(t)= xa(t) 
==
二、抽样定理
要想抽样后能不失真的还原出原信号,则抽样频率必须大于两倍的信号谱最高频率即，这就是抽样定理。
对连续信号进行等间隔采样形成采样信号，采样信号的频谱是原连续信号的频谱以采样频率为周期进行周期性的延拓形成的。
 =    (1.5.5)
三、抽样的恢复
如果满足抽样定理，则抽样后不会产生频谱混叠，故将通过如图所示的理想低通滤波器，就可得到信号频谱。
虽然理想低通滤波器是不可实现的，但在一定精度范围内可以用可实现的滤波器来逼近
下面讨论如何由抽样值来恢复原来的模拟信号。即通过H(jW)系统的响应特性。理想低通滤波器的冲激响应为
由与h(t)的卷积积分，即得理想低通滤波器的输出为
＝
＝
这就是内插值公式，即由信号的抽样值经此公式而得到连续信号，而称为内插函数，如图所示，在抽样点mT上，函数值为1，其余抽样点上，函数值为0。
在每个抽样点上，只有该点所对应的内插函数不为零，这使得各抽样点上信号值不变，而抽样点之间的信号则由各加权抽样函数波形的延伸叠加而成，如下图所示。这个公式说明了只要抽样频率高于两倍信号最高频率，则整个连续信号就可完全用它的抽样值来代表，而不会丢掉如何信息。这就是抽样定理的意义。
总结：如果序列是通过对模拟信号采样得到的，有关系：x(n)=xa(nT)，即序列值对于对模拟信号的采样值，或者说对于采样信号在t=nT时的幅度。
例：= sin(Wt)，理想抽样后，
x(n)= =sin(WnT)= sin(nω0)
∴ω0=WT数字域频率与模拟角频率之间的关系。
ω0=WT=W/fs=2pf/fs  2. 4时域离散信号的傅立叶变换与模拟信号傅立叶变换之间的关系
连续信号的傅立叶变换及反变换公式如下：
=
理想抽样后的抽样信号为，
= xa(t) p(t) =
则抽样信号的傅立叶变换
 = 
离散时间信号x(n)=xa(nT)，x(n)的傅立叶变换为：(2.2.1)
抽样信号的傅立叶变换与有什么关系？
可以证明，也可写成：
=，对照
 = ，都是以周期进行周期延拓。
画时，以w为横轴，以周期进行周期延拓。
画时，以W为横轴，以周期进行周期延拓。
坐标轴之间的对应关系如下图所示。
在一些文献中经常使用归一化频率，因为，和都是无量纲量，刻度是一样的。
所以数字频率0、2p处是低频，p附近代表高频。
当抽样频率是信号最高频率4倍时，最高频率所对应的数字频率为。

尚戈继：连续时间滤波器Matlab仿真实例​zhuanlan.zhihu.com
本文是这篇文章的后续，讨论了两类离散时间滤波器（无限冲激响应IIR和通过Kaiser窗函数法得到的有限冲激响应FIR）在Matlab中的实现。进行了对比分析。给出了代码
本篇目录：
待处理信号
离散时间滤波器
滤波器指标
IIR滤波器设计
FIR滤波器设计
FIR和IIR的效果对比
IIR和FIR的对比
完整代码
参考
待处理信号
与连续时间情况一样，待处理的信号是gamme.m函数生成的一段依次从Do、 Re、Mi、...一直到升Do的音频，每个频率持续一秒钟。
function [gam,t] = gamme(duree,fe)

% [do re mi fa sol la si doo]
freqnotes = [262 294 330 349 392 440 494 523];
t = 0:1/fe:duree;

gam = [];
for note=1:8, 
    gam = [gam, sin(2*pi*freqnotes(note)*t)];
end
N = length(gam);
t = (0:N-1)/fe;
调用函数，生成音频，并使用matlab的soudnsc函数试听
clear all
fe = 8192 ; %采样频率
duree = 1; %每个音符持续1秒
[sig,t] = gamme(duree,fe);
soundsc(sig, fe); % listen
对原始信号进行频率分析
fmin = 250;
fmax = 550;

f = [0:length(sig)-1]*(fe/length(sig));
S = fft(sig)/fe;

figure()
plot(f, abs(S));
axis([fmin fmax 0 0.5])
离散时间滤波器
滤波器指标
与连续时间情况不同，此时我们希望除去FA这个音，此时的滤波器指标
类型：陷波滤波器
通带容许的最小增益：-1dB
阻带增益：-40dB
过渡带中心频率：340Hz和360Hz(【Mi和Fa的中点】和【Fa和Sol的中点】)
过渡带宽度：10Hz 
IIR滤波器设计
0、IIR滤波器介绍
这里引用一段话[1]
历史上,随着数字信号处理领域的兴起,离散时间IIR滤波器的设计依赖于将连续时间滤波器变换成满足预定指标的离散时间滤波器的方法。该方法过去是,并且现在仍然是一种合理的方法, 其理由如下: 
•连续时间IIR滤波器的设计技巧十分成熟,并已取得许多有用的成果.因此可以方便地利用这些为连续时间滤波器推导出的设计方法。
•许多有用的连续时间滤波器设计方法有比较简单的完整设计公式。因此,以这种标准的连续时间IIR滤波器设计公式为基础的离散时间IIR滤波器的设计方法实现起来十分简单。
•当把完全适用于连续时间IIR滤波器的标准逼近方法直接用于离散时间IIR滤波器时,并不能得出简单的完整设计公式,因为离散时间滤波器的频率响应是周期的,而连续时间滤波器则不然。
1、常量定义
fcb = 340/(fe/2); %低频段的
fch = 360/(fe/2); %高频段的
deltaf = 10/(fe/2);
Wp = [fcb-deltaf/2,fch+deltaf/2]; %pass 通带的两个边界
Ws = [fcb+deltaf/2,fch-deltaf/2]; %stop 阻带的两个边界

Ap = 1;
Aa = 40;
这里Wp和Ws的定义见Cheby的文档 Chebyshev Type II filter order 
在离散的情况下频率需要做出对应的变换，即 离散频率=实际频率/(fe/2)，其中fe代表信号采样频率
2、计算滤波器阶数
这里我们使用的连续滤波器原型是Chebyshev Type II filter，在通带无抖动，在阻带有抖动
当然也可以选用其他类型的滤波器，方法类似，在连续时间滤波器Matlab仿真实例中介绍过了
[nd wnd]= cheb2ord(Wp,Ws,Ap,Aa) %这里不写‘s’就是离散的 %命名中的‘d’代表discrete
输出
注意这里得到的两个wnd是通带和过度带的交界位置频率。
3、计算滤波器参数
[numd, dend] = cheby2(nd,Aa,wnd,'stop')
可以画出这个滤波器的频域响应图
f = linspace(200,550,1000);
freqz(numd,dend,f,fe)
得到的滤波器很好地满足了设计指标，在阻带的增益小于-40（标的那条线），在增带的增益大于-1dB(图中看不太出来)
也可以看到IIR的缺点：相位是非线性的
滤波效果放在下下节，和FIR进行对比
FIR滤波器设计
0、FIR滤波器介绍
通常设计IIR滤波器采用的方法都是由连续时间系统到离散时间IIR系统的变换演变而来的。与之对应，离散FIR滤波器的设计方法以直接逼近所需离散时间系统的頻率响应或脉冲响应为基础。
设计FIR滤波器的最简单的方法称为窗函数法，核心思想在于：根据预期需求，在频域内设计一个“窗“函数--->得到该”窗“的脉冲响应，(非因果和无限长的)--->截断该理想脉冲响应，得到预期的FIR滤波器。
最容易想到的窗函数自然就是分段恒定的窗，如上图所示，但是，这种窗会引起吉布斯现象，即截断得到的FIR滤波器会随着有限响应长度（即截断的长度）增加，非一致地收敛到窗函数，而这种非一致收敛并不是我们想要的，具体来讲，这时滤波器并不能有效地滤掉不想要的频率成分，无论增加多少有限响应的长度。
因此，实际中并不会使用上图所示的窗，而是别的窗，如Bartlett、Hanning、Hamming、Blackman等，这些窗口的两端平滑地缩小到零，吉布斯现象会被减弱，但代价是主瓣加宽，过渡带的带宽增大，因此为了达到相同的性能指标，需要更高的阶数（即长度）。
这里使用的窗函数是Kaiser窗，与上图所示的窗不同，Kaiser窗有两个自由度，它是一种近似最佳的窗函数。
1、常量定义
这里不用做向离散的频率变换
fcb = 340;
fch = 360;
deltaf = 10;
Wp = [fcb-deltaf/2,fch+deltaf/2];
Ws = [fcb+deltaf/2,fch-deltaf/2];

Ap = 1;
Aa = 40;
2、滤波器系数计算
[n,wc,beta] = kaiserord([fcb-deltaf/2 fcb+deltaf/2 fch-deltaf/2 fch+deltaf/2]...
    ,[1 0 1],[1-10^(-Ap/20) 10^(-(Aa)/20) 1-10^(-Ap/20)],fe)
[n,Wn,beta,ftype] = kaiserord(f,a,dev,fs)
f：频率区间，a：不同区间的预期增益，[1 0 1]代表着[通 阻 通]
dev：不同区间内相较预期增益的容许偏差 我们的性能指标是以分贝的形式给的，这里要转换成小数 
fe：采样频率
3、滤波器
使用firl截取窗函数的脉冲响应，获得FIR滤波器
h = fir1(n,wc,'stop',kaiser(n+1,beta));
Kaiser窗函数的设置保证了得到的脉冲是因果并有广义线性相位特性的。因此滤波器的传递函数的分子系数num=h，分母系数den=1
同样的，画出该滤波器的频率响应
f = linspace(200,550,1000);
freqz(h,1,f,fe) %num = h den = 1
FIR和IIR的效果对比
频域响应上
figure()
hold on 
FIR = freqz(h,1,f,fe); %FIR
IIR = freqz(numd,dend,f,fe); %IIR

hold on
grid on 
plot(f,20*log10(abs(IIR)))
plot(f,20*log10(abs(FIR)))
legend('IIR','FIR')
Tips: FIR的误差大小约在3dB左右，因此，如果想让阻带的增益更加精确地小于-40dB，一个常用的技巧是在计算Kaiser窗函数系数时预先补偿这3dB，即
[n,wc,beta] = kaiserord([fcb-deltaf/2 fcb+deltaf/2 fch-deltaf/2 fch+deltaf/2]...
    ,[1 0 1],[1-10^(-Ap/20) 10^(-(Aa+3)/20) 1-10^(-Ap/20)],fe)；
% 原本是10^(-(Aa)/20)，改成10^(-(Aa+3)/20)
就能得到
滤波效果上
filter：调用离散时间滤波器处理信号
sigf = filter(numd,dend,sig); %IIR
sigf2 = filter(h,1,sig); %FIR %num = h den = 1
然后看输出的信号的频域成分
FFTplot函数的定义详见连续时间滤波器Matlab仿真实例文末，其实很简单，就是画出一个信号的傅里叶变换图，并放大到特定区间。
figure()
subplot(3,1,1)
FFTplot(sig, fe,200, 550) 
title('raw')

subplot(3,1,2)
FFTplot(sigf, fe,200, 550)
title('IIR')

subplot(3,1,3)
FFTplot(sigf2, fe,200, 550)
title('FIR')
 
差不多，但是输出信号在时域上则有不同
figure()
subplot(3,1,1)
plot(t, sig)
title('raw')
ylim([-2 2]);

subplot(3,1,2)
plot(t, sigf)
title('IIR')
ylim([-2 2]);

subplot(3,1,3)
plot(t, sigf2)
title('FIR')
ylim([-2 2]);
可以看到FIR有时延的效果。这源于它的广义非线性相位特性
IIR并没有时延的效果，IIR的相位时非线性的，得到的信号相位信息会被干扰，因此IIR适用于那些对相位信息不敏感的信号上。
IIR对信号相位信息的干扰，体现在
IIR处理得到的信号的”毛刺“更多一点（第4、5、6秒处）
IIR的通带和阻带的边界的滤波效果更差一点（第3秒处）
这种对相位信息的干扰，对于我们这里处理的相对比较简单的信号影响十分小。
可以试听一下滤波得到的声音
%不要三行命令同时运行，一次只运行一个
%如果电脑有杜比全景声（Dolby Atmos）的话，暂时关掉再播，否则听到的声音有失真
soundsc(sig, fe); % RAW listen
% soundsc(sigf, fe); % IIR listen
% soundsc(sigf2, fe); % FIR listen
IIR和FIR的对比
标准
FIR
IIR
相位
可控
不可控
复杂度
低
相对较高
稳定性
恒稳定
不恒稳定
阶数或参数个数
高
低
精准程度
普通
高
值得注意的是，虽然FIR在直观上是比IIR更加简单实现的，但是它的精确程度比IIR差，因此，对于同样的滤波器设计指标，FIR滤波器所要求的阶数可能比IIR滤波器高5-10倍，结果，成本较高，信号延时也较大；如果按线性相位要求来说，则IIR滤波器就必须加全通网络进行相位校正，同样要大大增加滤波器的阶数和复杂性。而FIR滤波器却可以得到严格的线性相位。[2]
IIR滤波器必须采用递归结构，极点位置必须在单位圆内，否则系统将不稳定。另外，在这种结构中，由于运算过程中对序列的舍入处理，这种有限字长效应有时会引入寄生振荡。相反，FIR滤波器主要采用非递归结构，不论在理论上还是在实际的有限精度运算中都不存在稳定性问题，运算误差也较小。此外，FIR滤波器可以采用快速付里叶变换算法，在相同阶数的条件下，运算速度可以快得多。[2]
以下三段来自[1]
选择IIR或是FIR滤波器取决于每种类型滤波器的优点在设计问题中的重要性。例如，IIR 滤波器具有可以用完整的设计公式来设计各种选频滤波器的优点。这就是说，一旦选定了选用哪种已知的逼近方法(即巴特沃思、切比雪夫或椭圆逼近)，则可以直接把技术指标代入一组设计方程来计算满足技术条件的滤波器的阶次，并得出数字滤波器的系数(或极点和零点)。这种IIR滤波器的非迭代计算程序。这些方法只限于设计选频滤波器，并只允许用于规定了幅度响应的场合。如果要得到其他形状的幅度响应，或需要逼近预定的相位响应或群延迟响应。则需要用算法设计法。
与此相反，FIR滤波器可以有精确的(广义)线性相位。但是对于FIR滤波器不存在完整的设计方程。虽然可以直接用窗函数数法，但为了满足预定的技术指标，有可能需要做一些迭代。与窗困数法相比。Parks-Mcclellan算法可以得出较低阶的滤波器，并且这两种方法的墟波器设计程序都很容易得到。而且。窗函数法和大多数算法设计法都有可能逼近较为任意的频率响应特性，但所遇到的困难要比在低通浦波器设计中遇到的稍大一些。。此外，FIR滤波器的设计问题要比IIR设计问题有更多的可控之处，因为对于FIR滤波器有适用于各种实际情况的最佳理论。
最后，在实现数字滤波器时还要考虑成本问题。通常将硬件的复杂性。芯片的面积或计算速度等作为衡量成本问题的因素。这些因素或多或少地直接与满足给定指标所需的滤波器阶次有关。在实际应用中，多相位实现的功效并没有表现出来，一般说来用IIR滤波器就能最有效地满足给定的幅度响应技术指标。但是在许多情况下，FIR滤波器的线性相位与它所带来的额外成本相比是非常值得的。
简单地说，FIR设计比较简单，IIR更加灵活，但也更加复杂。
完整代码
%信号生成
clear all
fe = 8192 ;
duree = 1;
[sig,t] = gamme(duree,fe);
soundsc(sig, fe); % listen
fmin = 250;
fmax = 550;

%常量定义
fcb = 340/(fe/2);
fch = 360/(fe/2);
deltaf = 10/(fe/2);
Wp = [fcb-deltaf/2,fch+deltaf/2];
Ws = [fcb+deltaf/2,fch-deltaf/2];

Ap = 1;
Aa = 40;


%%
%IIR 
%滤波器参数

[nd wnd]= cheb2ord(Wp,Ws,Ap,Aa)
[numd, dend] = cheby2(nd,Aa,wnd,'stop')

f = linspace(200,550,1000);
figure()
freqz(numd,dend,f,fe)
grid on
%滤波
sigf = filter(numd,dend,sig);

%%
%FIR 
%常量定义
fcb = 340;
fch = 360;
deltaf = 10;
Wp = [fcb-deltaf/2,fch+deltaf/2];
Ws = [fcb+deltaf/2,fch-deltaf/2];

Ap = 1;
Aa = 40;

%滤波器参数
E_comp = 3; % RIF最后在阻带会有3dB左右的误差，可以通过预先补偿的方式消除
[n,wc,beta] = kaiserord([fcb-deltaf/2 fcb+deltaf/2 fch-deltaf/2 fch+deltaf/2]...
    ,[1 0 1],[1-10^(-Ap/20) 10^(-(Aa+3)/20) 1-10^(-Ap/20)],fe)
% dev is a vector the same size as a that specifies the maximum allowable error...
%     or deviation between the frequency response of the output filter and its...
%     desired amplitude, for each band. The entries in dev specify the passband...
%      ripple and the stopband attenuation.
h = fir1(n,wc,'stop',kaiser(n+1,beta));
freqz(h,1,f,fe);

sigf2 = filter(h,1,sig);

%%
%画图

figure()
hold on
FIR = freqz(h,1,f,fe);
IIR = freqz(numd,dend,f,fe);

hold on
grid on
plot(f,20*log10(abs(IIR)))
plot(f,20*log10(abs(FIR)))
legend('IIR','FIR')

figure()
subplot(3,1,1)
FFTplot(sig, fe,200, 550)
title('raw')

subplot(3,1,2)
FFTplot(sigf, fe,200, 550)
title('IIR')

subplot(3,1,3)
FFTplot(sigf2, fe,200, 550)
title('FIR')

figure()
subplot(3,1,1)
plot(t, sig)
title('raw')
ylim([-2 2]);

subplot(3,1,2)
plot(t, sigf)
title('IIR')
ylim([-2 2]);

subplot(3,1,3)
plot(t, sigf2)
title('FIR')
ylim([-2 2]);

%%
%试听
% soundsc(sig, fe); % RAW listen
% soundsc(sigf, fe); % IIR listen
% soundsc(sigf2, fe); % FIR listen
参考
^ab离散时间数字信号处理，奥本海姆 https://book.douban.com/subject/26307919/
^abhttps://rf.eefocus.com/article/id-332067

【学习笔记】matlab进行数字信号处理（一）生成信号及信号的时域频域分析

【学习笔记】matlab进行数字信号处理（二）信号的相关分析及幅值分析

【学习笔记】matlab进行数字信号处理（三）数字滤波技术

【学习笔记】matlab进行数字信号处理（四）信号的时频域分析
b站视频地址：https://www.bilibili.com/video/BV18E411f7ZQ?p=16&spm_id_from=pageDriver
第二章 信号函数与标准信号
2.1 概述
信号发生器产生信号，可用信号函数生成正弦波、方波、三角波、锯齿波、白噪声、脉冲信号、阶跃信号、斜波信号、加速度信号
2.2信号函数和标准信号


连续的正弦波要进行离散化，才能进而的进行数字信号处理，将连续的时间t用n倍的时间间隔表示



改错：dt=1/11025











matlab中生成三角波和锯齿波用同一个函数sawtooth










2.3 信号发生器设计




b站视频P20：讲解如何生成信号发生器和电子琴（包含GUI设计）
第三章 信号的时域分析
3.1概述

3.2信号波形参数识别






画一条零线，可求出信号的周期和初相位

相邻两过零点时间差就是周期，第一个过零点位置与周期的比例可算出初相位

过零检测法求周期和相位很常用






3.3 信号的数字微分/积分







第四章 信号的频谱分析

4.1频谱分析的概念




可以从频谱图看出一个信号的频率构成，和不同频率不分的强弱（幅值）

从波形，可以看出信号随时间变化的幅值强弱

与波形相比，频谱可以更直观的

波形：幅值随时间的变化情况；频谱：幅值随频率的变化情况







与波形相比，频谱的抗干扰能力更强



频谱分析分解出的正弦波，往往能找到明确的物理意义，比如齿轮转动的频谱，每一个谱线都能对应到一个机械零件上去
4.2 周期信号的频谱分析（FFT代码和作业）




任意两个波形相乘，积分等于0说明两个函数正交



















功率谱是幅频谱的平方，反映了每个频率分量的能量大小



周期信号的频谱图：

谐波性是指频率成分一定是基频的倍数，如基频是f0，则谱线一定智能出现在f0的整数倍上

收敛性是指谐波次数越高，幅值越小





fft使用条件：数据长度必须是2的次方



合成信号——FFT——实频谱和虚频谱——幅值谱和相位谱



【视频P27，11分位置】注意：使用FFT之后，并不需要显示负频谱的部分，因为正频谱部分就可以显示出有用的信息，负频谱部分的信息是冗余的



有时功率谱上某些频谱分量的幅值显示的不明显，小能量信号会被大能量信号掩盖，所以通常采用对数功率谱，能将小的信号放大，大的信号压缩，能在一张图上表现出不同大小的频率成分




4.3 数字信号的频谱计算方法



截取后再延拓的信号会在连接点处产生跳变，在波形上显示的是产生跳变，在频谱上产生的是出现能量泄露












4.4 FFT中的能量泄露和栅栏效应


快速傅立叶变换和离散傅立叶变换的计算结果是相同的，但是速度快很多





栅栏效应生成的原因：为了提高计算效率，在使用FFT计算频谱时，在0到二分之一的采样频率区间内，不是每个频率点都计算，而是按照一定的频率间隔抽样计算，频率间隔等于采样频率除以数据长度

在取样的位置和谱峰的位置不重叠时，便会带来一个谱峰误差，这种误差引起的反应，叫做栅栏效应。包含两个误差，一个是幅值高度的误差，一个是最大频率的误差。



谱峰越尖锐，产生误差的可能性就越大





由于截断带来的能量泄露，这时会减小栅栏效应的误差，可以通过控制截断函数的形状，来调整能量泄露形状，

左侧图为矩形窗能量泄露函数，以中心频率为主的成份叫做主瓣，比较窄比较尖，越尖锐造成的栅栏效应的误差就越大，还有其他很高的旁瓣，但理想情况下只希望读出主瓣处的能量泄露，但这种理想的窗函数并不存在，实际上只能去逼近它。


直接截断，然后周期延拓，如果截断的不是整周期的位置，会产生一个跳变，此时加一个汉明窗，用这样的窗函数和原来的信号相乘，这样会将原信号在接头处的幅值被压缩成0，以这种方式进行周期延拓，无论原来的信号怎样，在接头处都不存在跳变（从波形上看）

（从频谱上看）相当于增宽主瓣，压缩旁瓣，因此可以通过这种方式抑制旁瓣的能量泄露，同时提高主瓣的宽度，使栅栏效应的误差减小







汉宁窗和矩形窗相比，主瓣变的更宽，而旁瓣变得更低



这是汉宁窗对波形进行截断的情况，截断后要进行修正，窗函数和原信号相乘再进行截断延拓，幅值修正后幅值增大了一倍，这样幅值增大后可以弥补两端衰减造成的能量损失







和矩形窗相比，主瓣宽一些，旁瓣窄一些



主瓣更宽一些，旁瓣更窄一些





平顶窗，主瓣相当宽





能量修正z=2*w1.*x，然后再进行FFT

加窗能克服能量泄露和栅栏效应带来的误差


4.5 非周期信号的频谱分析


工程上大量的信号都是非周期信号，在数字领域，由于截断，不管原信号是不是周期信号，都会转换成周期信号



1.如果数据长度不够长，可以通过补0的方式提高频谱的频率分辨率

eg：假如原来数据长度是1024，再补1024个0，频谱分析时，精度会高一倍

2.细化：如果想看清某一部分的频谱，使用ZOOM-FFT可以将观测的视角集中在一个频率

3.主要针对受噪声干扰频谱，计算时谱当中有很多的随机干扰，此时将频谱多次累加再除以N，这样干扰也会减小到N分之一

4.用于修正栅栏效应误差

5.传感器输出的都是实信号，但是matlab里的FFT输入的是复信号，因此针对虚部的计算浪费了，采用实信号FFT计算技术，可以将计算速度提高一倍
4.6频谱分析应用




结果会发现，信号波形的幅值是1，倒是频谱的幅值不是1，产生差异的原因就是能量泄露和栅栏效应
大作业：声音信号采集和频谱分析程序设计
首先通过计算机上的脉冲，采集一段声音信号，然后进行FFT变换，再画出信号频谱

1.新建GUI文件

2.安放控件，修改控件属性

3.使用定时器，实现连续采样