第3章 小波分析基础

众所周知，傅立叶变换一直是时频分析的主要方法，并且广泛地应用于工程实践。傅立叶变换将信号只在频域上展开，因为信号频率的信息对其是非常重要，所以对于某些应用来说是很恰当的。
同时，傅立叶变换又是小波变换的基础，因此，有必要先了解傅立叶变换的理论及其应用，故本书将此章作为学习小波分析的起点。学习好本章内容的同时，读者也将会了解到傅立叶变换等一些基本时频变换方法在MATLAB中的应用以及信号处理的方法。
本章的主要安排是：前两节主要介绍傅立叶变换及其在MATLAB中的应用，主要包括离散傅立叶变换DFT和快速傅立叶变换FFT等；另外，介绍了Z变换及其在MATLAB中的应用；最后，针对一些滤波器例举了一些常见的实例。
学习目标：
(1)了解傅立叶变换的基本概念
(2)熟练掌握一二维傅立叶变换的函数使用
(3)熟练掌握函数residuez
(4)熟练掌握滤波器
3.1 一维傅立叶变换及其应用
传统的信号分析是建立在傅立叶(Fourier)变换的基础上的，并且在具体工程实践中，研究者对其进行了推广甚至根本性的革命，如Z变换、短时傅立叶变换、小波变换等。
3.1.1 一维傅立叶变换
傅立叶变换是频谱理论的基础。较早时，往往利用傅立叶级数即调和分析的方法来分析频谱特征。随着计算机的发展，已经逐渐开发出这些方法的程序。
从实用的角度而言，人们考虑傅立叶分析通常是指傅立叶变换和傅立叶级数。对一维信号 f (t)来说，其连续地傅立叶变换定义为：
对于参数ω的任何实数值都存在，则称F(ω)为f (t)的傅立叶变换，或称傅立叶积分。函数F(ω)一般为复数，可以写成：
其中，A(ω)称为f (t)的傅立叶谱(振幅谱)，Φ(ω)为相位谱，A2
 (ω)为能量谱。
而它的傅立叶逆变换定义为：
为了计算傅立叶变换，需要用数值积分，即取 f (t)在时间轴上的离散点上的值来计算这个积分。接下来将进一步作具体介绍。
3.1.2 一维离散傅立叶级数
若t是定义在时间轴上的离散点，则x(t)为离散时间信号，表示为x(nT)，其中，T表示相邻两个点之间的时间间隔，也称为抽样周期，n 取整数并代表时间的离散时刻。一般可以把T归一化为1，这样x(nT)可简记为x(n)，又称为离散时间序列。
系统是由若干相互作用和相互依赖的事物组合而成的并具有特定功能的整体。系统分析的着眼点是分析系统的输入和输出的关系，而不涉及系统的内部情况。
一个离散系统可以抽象为一种变换或者一种映射，即y(n)=T[x(n)]，其中，T代表变换。在数字信号处理中，所研究的系统基本上都是线性时不变系统(LSI系统)，其输入/输出关系可以通过冲激响应h(n)表示：
其中，x(n)和h(n)作卷积运算，MATLAB提供了求卷积函数conv。有关线性时不变系统中涉及的重要概念可参看信号与系统的经典教材。
在MATLAB中，可用一个向量来表示一个有限长度的序列(由于内存的限制，不可能表示一个任意无限序列)，示例代码设置如下：
x(n)=[1，3，2，-3，6，7，0，-8，9]
MATLAB对下标的约定从1开始递增(这里与C语言有区别)，如果要包含采样时刻的信息，则需要两个向量来表示。
为了分析的方便，数字信号处理理论中定义了一些基本的典型的序列，其定义和MATLAB表达式如下。
1.单位取样序列
该表达式可利用zeros函数来实现，代码设置：
x=zeros(1，N);
x(1)=1;
2.单位阶跃序列
该表达式可利用ones函数来实现，代码设置如下：
x=ones(1，N)
3.实指数序列
该表达式通过MATLAB实现，代码设置如下：
n=0:N-1;
x=a.^n;
4.复指数序列
该表达式通过MATLAB实现，代码设置如下：
n=0:N-1;
x=exp(a+j*w0)*n);
5.正(余)弦序列
该表达式通过MATLAB实现，代码设置如下：
n=0:N-1;
x=cos(w0*n+q);
6.随机序列
在MATLAB中提供了两类(伪)随机信号：rand(1，N)产生[0，1]上均匀分布的随机矢量；randn(1，N)产生均值为0、方差为1的高斯随机序列。其他分布的随机数可通过上述随机数的变换而产生。
7.周期序列
如果信号在频域上是离散的，则该信号在时域上就是周期性的函数；反之，在时域上离散的信号在频域上必然表现为周期性的频率函数。
可以得出一个一般规律：一个域的离散必然造成另一个域的周期延拓。这种离散变换本质上都是周期的。
对于周期为N的离散时间信号序列
 ，其中k为任意整数，由于在Z平面上没有任何收敛区域，所以不能进行Z变换，但是可以用傅立叶级数来表达，其基波频率为
 ，用复指数表示为
 ，第k次谐波为

 ，所以有
 ，可得离散傅立叶级数公式如下：
在上式中求和号前所乘的系数1/N是习惯上采用的常数，
 是k次谐波的系数：
习惯上也常采用符号
 ，这样离散傅立叶级数对可以表示为：
3.1.3 一维离散傅立叶变换及应用
周期序列实际上只有有限个序列值有意义，因此它的许多特性可以沿用到有限长序列上，对于一个长度为N的有限长序列x(n)，以x(n)为主值序列并以N为周期进行延拓得到周期序列
 ，即：

 ；n为其余值时，x(n)为0。
由离散傅立叶级数公式可以得到有限长序列x(n)的离散傅立叶变换公式：
下面介绍一个离散傅立叶变换的实例。
【例3.1】在MATLAB中 离散傅立叶的实现。
保存在eg3_1.m中，设置MATLAB程序代码如下：
N=16;　　　　%序列长度
n=0:N-1;　　　　%时域取样
xn=sin(pi*n/4);　　　%产生序列
k=0:N-1;　　　　%频域取样
wn=exp(-j*2*pi/N);
nk=n'*k;
wnnk=wn.^nk;
xk=xn*wnnk　　　%计算DFT
%画图figure(1)
stem(n，xn)
figure(2)
stem(k，abs(xk))
结果如下，有限序列的DFT和序列的DFT分别如图3-1和图3-2所示。
xk=
Columns 1 through 5
4.0137　　　8.5406 -7.6853i　-4.1658+ 4.5586i　-3.0693+ 3.2036i　-0.8112+ 0.7734i
Columns 6 through 10
-0.7071 + 0.5412i -0.6552 + 0.3857i -0.6266 + 0.2676i -0.6104 + 0.1696i -0.6021 + 0.0824i
Columns 11 through 15
-0.5995 - 0.0000i -0.6021 - 0.0824i -0.6104 - 0.1696i -0.6266 - 0.2676i -0.6552 - 0.3857i
Columns 16 through 20
-0.7071 - 0.5412i -0.8112 - 0.7734i -3.0693 - 3.2036i -4.1658 - 4.5586i 8.5406 + 7.6853i图3-1 有限序列的DFT
图3-2 序列的DFT
在进行快速傅立叶变换FFT的操作时可以调用内部函数fft，这样速度比较快。DFT和FFT在信号处理中有着重要的应用，可以进行卷积运算，实现线性时不变系统等。
说明：此算例保存为eg3_1.m，无特殊说明以后章节的算例在MATLAB中的保存均以类似命名方式命名。
3.1.4 一维快速傅立叶变换及应用
快速傅立叶变换是离散傅立叶变换的快速算法，其应用极大地推动了DSP理论和技术的发展。下面介绍FFT的基本思想和实现。
FFT的基本思想是在N点序列的DFT为：
由于系数
 是一个周期函数：
 ，且是对称的：
FFT正是基于这样的基本思想而发展起来的。下面介绍一个快速傅立叶变换的实例。
【例3.2】计算正弦信号波的快速傅立叶变换，并画出频率与相位图。
结果保存在eg3_2.m中，代码设置如下：
t=(0:99)/100;　　　　　% 时间向量
x=sin(2*pi*15*t)+ sin(2*pi*40*t);　　% 正弦信号波
y=fft(x);　　　　　%调用快速傅立叶变换函数
p=unwrap(angle(y));　　　　%相位
f=(0:length(y)-1)'/length(y)*100;　　%频率向量
plot(f，p)
结果如图3-3所示。
图3-3 FFT相位图
【例3.3】设x(n)长度为N1
 ，h(n)长度为N2
 ，计算两者的线性卷积，其步骤如下：
(1)将x(n)和h(n)通过增加零值延长到N=N1
 + N2
 −1 ；
(2) 利用FFT计算x(n)和h(n)各自N 点的DFT；
(3)计算Y(k)=X(k)H(k)；
(4) 利用IFFT 计算Y(k)便可得到卷积结果y(n)。
MATLAB程序代码如下：
x=[5，2，3，8，9，1];
h=[8，2，4，1，5，8，2，3，9，5];
lenx=length(x);
lenh=length(h);
N=lenx+lenh-1;
xk=fft(x，N);
hk=fft(h，N);
yk=xk.*hk;
y=real(ifft(yk));
plot(y);
xlabel('n');
ylabel('y(n)');
title('x(n)*h(n)');
grid
运行结果如图3-4所示。
图3-4 应用FFT 计算线性卷积
3.2 二维傅立叶变换及其应用
上节介绍了一维傅立叶变换及其应用，本节将介绍二维傅立叶变换及其应用。二维离散傅立叶变换的定义如下：
对于二维离散傅立叶变换，MATLAB提供了一个函数，其调用格式如下：
Y=fft2(X)
其中：X是二维信号；Y是返回与X同维的数组。
下面将应用此函数，利用二维离散傅立叶变换在MATLAB中实现图像变换技术。
【例3.4】对矩阵A进行二维傅立叶变换。MATLAB代码设置如下：
load imdemos saturn;
subplot(1，2，1);　　　　%显示图像
imshow(saturn);
%进行傅立叶变换
B=fftshift(fft2(saturn));
subplot(1，2，2);　　　　%显示变换后的系数分布
imshow(log(abs(B))，[])
colormap(jet(64))
colorbar;
结果如图3-5所示。
图3-5 二维傅立叶变换图
接下来介绍一个二维傅立叶逆变换在MATLAB应用的例子。
【例3.5】在图像text.png中定位字符“a”。
傅立叶变换还能够用来分析两幅图像的相关性，相关性可以用来确定一幅图像的特征，在这个意义下，相关性通常被称为模板匹配。
将包含字母“a”的图像与text.png图像进行相关运算，也就是首先将字母a和图像text. png 进行傅立叶变换，然后利用快速卷积计算字母 a和图像 text.png的卷积(其结果如图3-6b所示)，提取卷积运算的峰值，如图3-6d所示的白色亮点，即得到在图像text.png中对字母“a”定位的结果。
MATLAB代码设置如下：
bw=imread('text.png');
%从图像中抽取字母a的图像
a=bw(32:45，88:98);
subplot(2，2，1)，
imshow(bw);
subplot(2，2，2)，
imshow(a);
C=real(ifft2(fft2(bw).*fft2(rot90(a，2)，256，256)));
subplot(2，2，3)
imshow(C，[])
%选择一个略小于C中最大值的值做为阈值
max(C(:));
thresh=60;
subplot(2，2，4)，
%显示像素值超过阈值的点
imshow(C>thresh)
结果如图3-6所示。
图3-6 图像特征识别示例

3.3 Z变换及其应用
Z变换是离散系统和离散信号分析与综合的重要工具，也是傅立叶变换的推广，在MATLAB中也有相应的函数。首先，了解Z变换的定义。
给定离散信号x(n)，其Z变换的定义为：
其中，z为复变量。如果 n的取值范围为−∞～+∞，则上式定义的Z变换称为双边Z变换；如果n的取值范围为0～+∞，则上式称为单边Z变换。
实际的物理系统的抽样响应h(n)在n<0时恒为0，因此对应的都是单边Z变换。
X(z)存在的z的集合称为收敛域(ROC)，即满足下式：
ROC由|z|决定，因此一般为环形。序列在单位圆上的Z变换就是序列的傅立叶变换，若x(n)=h(n)，则该变换就是系统的频率响应。
Z变换的性质如下。
(1)线性性质。
若Z[x1
 (n)]=X1
 (z)，ROC:R1
 ；Z[x2
 (n)]=X2
 (z)，ROC:R2
 ，则有：

 ，ROC为R1
 和R2
 的公共部分。
(2)时移性质。
(3)频移性质。
(4)折叠。

 ，ROC:ROC与ROCx
 相反，即ROCx
 的逆。
(5)复共轭。
Z(x*(n))=X*(z*)，ROC:ROCx
(6)Z域微分。

(7)序列相乘。

 的逆，其中c为公共ROC中包围原点的闭合围线。
(8)序列卷积。
Z变换提供了任意序列在频域的表示方法，但它是连续变量 z的函数，因此无法直接利用计算机进行数值计算。为了使用MATLAB，必须截断序列，得到有限个点的表达式，这就产生了离散傅立叶级数(DFS)、离散傅立叶变换(DFT)和计算量小的快速傅立叶变换(FFT)。
对于Z变换为X(z)的序列，MATLAB的表示是通过X(z)的系数实现的，Z变换经常用到deconv、residuez和freqz等函数。
函数residuez的调用格式如下：
[r，p，k]=residuez(b，a)
其中，b和a分别是多项式的系数向量。
【例3.6】计算分数表达示的Z变换。
MATLAB程序代码如下：
>>b=1;
>>a=poly([1，-0.6]);
>>[R，P，C]=residuez(b，a)
结果代码为：
r=
0.6250
0.3750
p=
3.0000
-0.6000
k=
[]
因此得到上述分数表达式的Z变换为：
信号处理工具箱里包含用于对上述理论实现的函数，其功能参见表3-1。
表3-1 常用变换函数归纳表
续表
3.4 滤波器
时间域中信号处理的重要方法之一是滤波，即将一个混合信号中的某一分量分享出来或是把它去掉。
3.4.1 连续滤波器
令滤波器的脉冲响应为h(t)，其傅立叶变换为H (ω)，则滤波器在时间域里的计算公式为：
其中：f (τ)为输入，g(t)为输出。
在频域中的计算公式为：
式中：G(ω)为g(t)的傅立叶变换，F(ω)为f (τ)的傅立叶变换。
任意信号通过滤波器后，其振幅和相位都要发生变形，这就是滤波作用。例如，输入一个正弦波：
f (t)的傅立叶变换为：