门爱东老师讲义第章PPT课件.ppt

1 FFT FastFourierTransform1965年 JamesW Cooley和JohnW Tukey在 计算数学 MathematicsofComputation 上发表了 一种用机器计算复序列傅立叶级数的算法 AnalgorithmforthemachinecalculationofcomplexFourierseries 论文 自此之后 新的算法不断涌现 一种是对N等于2的整数次幂的算法 如基2算法 基4算法 另一种是N不等于2的整数次幂的算法 例如Winagrad算法 素因子算法 3 4FFT 快速离散傅里叶变换 Dr JamesW Cooley WorkedatIBMWatsonResearchCenterinYorktownHeights N Y AfterhisretirementfromIBMin1991 hejoinedtheElectricalEngineeringDepartmentattheUniversityofRhodeIsland JohnW Tukey 1915 2000 2 3 4FFT 直接计算DFT的运算量分析 N点有限长序列x n 的DFT变换对的定义为 其中 假设x n 是复序列 同时X k 一般也是复数 3 3 4FFT 直接计算DFT的运算量分析 复乘的加法次数 复加的加法次数 4 如N 512 1024和8192时 DFT的乘法运算N2 5122 218 262144 26万次 N2 10242 220 1048576 105万次 N2 81922 226 67108864 6千7百万次 对于大N 在实际中是不能接受的 无法 实时 应用DFT 一般地 在计算机中 一次加法比一次乘法所需时间要短 在DSP中 由于乘法用特殊的硬件电路专门完成 因此乘法和加法所需机器周期相同 Cooley与Turkey提出的FFT算法 大大减少了计算次数 如N 512时 FFT的乘法次数约为2000次 提高了约128倍 而且简化随N的增加而巨增 因而 用数值方法计算频谱得到实际应用 3 4FFT 直接计算DFT的运算量分析 5 3 4FFT WNkn的性质 6 以4点DFT为例 直接计算需要 42 16次复数乘 而按周期性及对称性 可以将DFT表示为 只需要1次复数乘 3 4FFT WNkn的性质 7 FFT算法分类 时间抽选法DIT Decimation In Time频率抽选法DIF Decimation In Frequency 3 4FFT 基本思想 基本思路 虽然存在不同的FFT方法 但其核心思想大致相同 即通过迭代 反复利用低点数的DFT完成高点数的DFT计算 以此达到降低运算量的目的 迭代 利用WNkn的周期性 特殊点和对称性 合并DFT计算中很多重复的计算 达到降低运算量的目的 低点数 将傅氏变换DFT分解成相继小的DFT计算 即N变小 而计算量与N2成正比 8 算法原理设序列点数N 2M M为整数 若不满足 则补零 将序列x n 按n的奇偶分成两组 N为2的整数幂的FFT算法称基 2FFT算法 即一组由偶数序号组成 另一组由奇数序号组成 3 4 1FFT 基2时间抽选法 算法原理 注意其长度为N 2 9 3 4 1FFT 基2时间抽选法 算法原理 10 偶数取样点DFT为 3 4 1FFT 基2时间抽选法 算法原理 奇数取样点DFT为 k的整个范围为0 N 1 而X1 k X2 k 是由N 2个样点形成的DFT x 2r 和x 2r 1 的长度为N 2 由这两个偶数和奇数N 2个时域样值可以计算出前N 2个DFT系数 也可以计算出后N 2个DFT系数 问题 这前后N 2个DFT有无关系 k在N 2 N 1 时 X1 k X2 k WN情况如何 11 3 4 1FFT 基2时间抽选法 算法原理 观察 则 12 因此 整个X k 的计算 可以分解为前 后半部分的运算 而只要求出前一半 就可以由上式求出整个序列 3 4 1FFT 基2时间抽选法 算法原理 同理 而 因此 13 上式表示为信号流程图 此信号流程图也称为蝶形流程图 3 4 1FFT 基2时间抽选法 算法原理 MorphoButterfly大闪蝶 南美洲 14 以N 8为例 其信号流程图 3 4 1FFT 基2时间抽选法 算法原理 偶数序列 奇数序列 15 3 4 1FFT 基2时间抽选法 算法原理 16 3 4 1FFT 基2时间抽选法 算法原理 同理 17 以N 8为例 其信号流程图为 3 4 1FFT 基2时间抽选法 算法原理 18 由于N 2M 这样逐级分解 直到2点DFT当N 8时 即分解到X3 k X4 k X5 k X6 k k 0 1 3 4 1FFT 基2时间抽选法 算法原理 19 每一次分解都是按时间域输入序列的奇偶性次序 分解为两个半长的子序列 所以称为按时间抽取法 仍以N 8为例 3 4 1FFT 基2时间抽选法 算法原理 m 0级 m 1级 m 2级 20 2m点DFT的基2时间抽选计算过程 可见 为什么引入FFT 出发点 考虑W因子的特点和性质 简化算法 DIT FFT算法 时间域输入信号逐级分解为奇偶序列 3 4 1FFT 基2时间抽选法 算法原理 上式中 21 例3 25使用基2时间抽选法FFT流图 计算下列数据的8点DFT x 4 3 2 0 1 2 3 1 3 4 1FFT 基2时间抽选法 算法原理 22 蝶形运算在FFT信号流图中每一级里节点都是成对存在的 计算时总是下节点的值乘以WNr 然后与上节点的值相加减 形成下一列两个节点值 这种计算的基本关系是 式中p q是上下节点的序号 从0开始 p q 0 1 N 2 每一级中每对节点称为对偶节点 对偶节点的间距为 注意 只有下节点乘以加权因子WNr 3 4 1FFT 基2时间抽选法 蝶形运算 23 同址计算 同位特性 3 4 1FFT 基2时间抽选法 同址计算 1 不同级数的节点序号不变从蝶形运算可以看出第m列的xm p xm q 经蝶形运算后 在第 m 1 列中的节点序号不变 即xm 1 p 中的p与xm 1 q 中的q仍分别是xm p xm q 中的p和q值 24 2 蝶形运算是自成独立单元蝶形运算自成独立单元 即xm 1 p xm 1 q 只与xm p xm q 有关 而与其他节点的值无关 同时xm p xm q 也不参与另外的蝶形运算 因此 就可把计算结果xm 1 p xm 1 q 放入计算前xm p xm q 的存储单元中 而不用再建新的存储单元 同址计算 同址计算所需存储量仅等于给定数据所需的存储量 这可大大节省存储单元 3 4 1FFT 基2时间抽选法 同址计算 25 输入倒序重排N点DFT分解为两个N 2点DFT 输入序列按奇偶分组 再分解 再奇偶重排 直到2点DFT 即FFT输出自然序列 输入序列x n 倒序重排 3 4 1FFT 基2时间抽选法 倒序重排 26 说明 首先把x n 中的n表示为二进制 假定N 8 则x n x n2n1n0 ni 0 1FFT的时间抽选法按n的奇偶分离而形成倒置重排的原理如上图所示 由此可以看出 时间抽选法FFT的输入倒序重排 输出结果为自然顺序 实际操作输入序列重排实际上就是完成二进制前后位序的相互交换 3 4 1FFT 基2时间抽选法 倒序重排 27 当N 2M时 共有M log2N级蝶形 每级N 2个蝶形 每个蝶形有1次复数乘法 2次复数加法 复数乘法 复数加法 比较DFT FFT 3 4 1FFT 基2时间抽选法 运算量 28 FFT算法上面讨论的FFT算法假定N为2的整数次幂 即N 2M 是基2FFT算法 当N Rv时 这样的算法叫做基RFFT算法 而当N R1v1R2v2R3v3时 叫做混合基算法 当N是一个高度合数时 可得到最有效的算法 最受欢迎也最易编程的算法是基2FFT算法 对于不能分解的质数 或者当N不是2的整数次幂时 可以在信号的末尾补0 使其成为高度合数或2的整数次幂 MatlabFFT实现Matlab提供了内建的X fft x N 函数来计算矢量x的DFT fft函数是机器码写成的 而不是以Matlab指令写成的 即不存在 m文件 因此它的执行速度很快 采用混合算法 若N为2的幂 则得到高速的基2FFT算法 若N不是2的幂 则将N分解成质数 得到较慢的混合基FFT 若N为质数 则fft函数采用原始DFT算法 3 4 1FFT 基2时间抽选法 Matlab实现 29 如果输入序列x的长度小于N 则在其填零使其成为N点序列 如果省略变量N 则DFT的长度即为x的长度 如果x为一个矩阵 则计算x中每列的N点DFT IDFT由ifft函数完成 它的特征与fft函数相同 3 4 1FFT 基2时间抽选法 Matlab实现 例3 26研究当1 n 2048时 fft函数的执行时间 这将展示不同的N的情形下 划分 组合的策略 解 atlab提供了两个函数来确定执行时间 clock函数读取瞬时时钟 etime t1 t2 函数计算时刻t1 t2之间所经历的时间 为了确定执行时间 产生长度为1至2048的随机矢量 计算它们的FFT 将计算时间存在一个数组里 最后画出执行时间相对于N的图 30 ComputationalComplexityofFFTusingMATLABNmax 2048 fft time zeros 1 Nmax forn 1 1 Nmaxx rand 1 n Uniformlydistributedrandomnumbers t clock fft x fft time n etime clock t 花费的时间endn 1 1 Nmax top max fft time plot n fft time axis 0 2500 0 top xlabel N ylabel TimeinSec title FFTexecutiontimes text 2100 top o N N text 2100 top 3 4 o N N 3 4 text 2100 top 2 o N N 2 text 2100 top 3 o N N 3 text 2100 top 4 o N N 4 text 2100 min fft time o N logN 3 4 1FFT 基2时间抽选法 Matlab实现 31 从图中可以看到 图中曲线分成几组 表示了N的不同组合情况 当N高度可合时 划分 组合策略是很有效的 例如N 2048时 执行时间约为0秒 N 2047时为0 05秒 N 2029时为0 22秒 在实际中必须谨慎使用高度合数 除特定需求外 实际中的一般选择是N 2M 当信号的取样点数不是2的整数次幂时 可以在信号的末尾补0 Zeropadding 外加的0不会影响信号的DTFT或DFT特性 当也要注意 填0后可能带来的 真实频率 模糊问题 见前面的例题 3 4 1FFT 基2时间抽选法 Matlab实现 32 DFT定义表示为矩阵形式为 3 4 2FFT 基2时间抽选法 矩阵表示 33 令 矩阵WN为 DFT矩阵简单地表示为 3 4 2FFT 基2时间抽选法 矩阵表示 34 FFT基2时间抽选法信号流图 N 8 1 W82 1 8点 2 4点 4 2点 x 0 x 1 x 2 x 3 x 4 x 5 x 6 x 7 x 0 x 4 x 2 x 6 x 1 x 5 x 3 x 7 X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7 1 1 1 1 1 1 1 1 1 1 1 W80 W80 W82 W80 W81 W82 W83 C 0 C 1 D 0 D 1 E 0 E 1 F 0 F 1 A 0 A 1 A 2 A 3 B 0 B 1 B 2 B 3 m 0级 m 1级 m 2级 倒置重排 3 4 2FFT 基2时间抽选法 矩阵表示 35 FFT矩阵形式表示为 1 输入矢量x与P8相乘 则只是输入的倒置重排 没有乘法操作 3 4 2FFT 基2时间抽选法 矩阵表示 36 2 第0级运算 2点DFT 3 4 2FFT 基2时间抽选法 矩阵表示 37 3 第1级运算 4点DFT 3 4 2FFT 基2时间抽选法 矩阵表示 38 4 第2级运算 8点DFT 3 4 2FFT 基2时间抽选法 矩阵表示 39 FFT算法看成是DFT矩阵W8的分解 这个因式分解对应于快速算法 因为矩阵F8 8 F8 4 F8 2 和P8的大部分元素为0 是稀疏矩阵 因此上述每个矩阵的乘法运算最多只需要8次复乘运算 而P8只是置换操作 没有乘法操作 不同的FFT算法对应于将DFT矩阵WN分解成不同的稀疏矩阵 3 4 2FFT 基2时间抽选法 矩阵表示 40 输入倒位序 输出自然序输入自然序 输出倒位序输入输出均自然序各级具有相同几何形状输入倒位序 输出自然序输入自然序 输出倒位序 3 4 2FFT 基2时间抽选法 其它形式的流图 41 3 4 2FFT 基2时间抽选法 其它形式的流图 42 3 4 2FFT 基2时间抽选法 其它形式的流图 43 3 4 2FFT 基2时间抽选法 其它形式的流图 44 3 4 2FFT 基2时间抽选法 其它形式的流图 45 3 4 2FFT 基2时间抽选法 其它形式的流图 46 算法原理 根据时间 频率的对偶性时间抽选法 是把输入x n 按奇偶分解成两个子序列 即N点x n 序列 N 2点子序列 而输出X k 是按自然顺序排列的 频率抽选法 是把输入x n 按照前后对半分开 而不是奇偶数分开 而输出X k 逐项分解成偶数点子序列和奇数点子序列 DFT变换表达式为 3 4 3FFT 基2频率抽选法 DIF FFT 如果将输入x n 按前后等分 即将求和分成两部分 范围分别为 47 3 4 3FFT 基2频率抽选法 算法原理 48 按k的奇偶将X k 分成两部分 3 4 3FFT 基2频率抽选法 算法原理 49 令 则X 2r 和X 2r 1 分别是x1 n 和x2 n 的N 2点DFT 记为X1 k 和X2 k 3 4 3FFT 基2频率抽选法 算法原理 50 3 4 3FFT 基2频率抽选法 算法原理 51 FFT基2频率抽选法信号流图 N 8 逐级分解 直到2点DFT 3 4 3FFT 基2频率抽选法 算法原理 52 基本蝶形运算不同DIT 先复乘后加减 W因子在上下节点都有体现 DIF 先减后复乘 W因子仅体现在下节点 3 4 3FFT 基2时间和频率抽选法的异同 53 运算量相同 同址计算 3 4 3FFT 基2时间和频率抽选法的异同 DIF输出重排 DIT输入重排 54 DIT和DIF的基本蝶形互为信号流图转置 DIT DIF 3 4 3FFT 基2时间和频率抽选法的异同 55 DIT DIFFFT算法互为转置 转置定理 1 1 8点 2 4点 4 2点 X 0 X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 0 X 4 X 2 X 6 X 1 X 5 X 3 X 7 x 0 x 1 x 2 x 3 x 4 x 5 x 6 x 7 1 1 1 1 1 1 1 1 1 W80 W82 W80 W81 W82 W83 C 0 C 1 D 0 D 1 E 0 E 1 F 0 F 1 A 0 A 1 A 2 A 3 B 0 B 1 B 2 B 3 m 0级 m 1级 m 2级 倒置重排 1 1 W80 W82 3 4 3FFT 基2时间和频率抽选法的异同 56 前面讲的都是基2的FFT算法 除此之外 还有基4的 基8的快速算法 原理和基2的类似 分解为4个交错的集合 相比基2的 可以进一步节约复数乘的次数 但是基8的和基4的差别不远 算法原理省略 自学 3 4 4FFT 基4时间抽选法 57 DFT和IDFT的定义 3 4 5FFT IDFT快速算法 IFFT DFT和IDFT的区别 因子W的指数相差一个负号 相差一个因子1 N FFT算法中的分组方式 排序方式以及蝶形运算结构都可用于IFFT算法的设计 而这就是可依据现有的FFT算法直接得出IFFT算法的原因 58 FFT和IFFT基本蝶形运算之间的关系设有序列x n 其DFT为X k 则IDFT为 在FFT的时间抽选法中 3 4 5FFT IDFT快速算法 IFFT 对于IFFT算法 输入是X k 和X k N 2 输出是X1 k 和X2 k 解上式可以得到X1 K X2 K 59 8点DIT IFFT算法 说明 1 分组过程是按时间序列x n 的奇偶性在时域上展开的 故称此法为时间抽选算法DIT IFFT 2 1 N的分解 N 2m 3 4 5FFT IDFT快速算法 IFFT 60 8点DIF IFFT算法 3 4 5FFT IDFT快速算法 IFFT 61 用FFT程序求IFFT的方法直接法 利用DIT DIF的FFT程序 改变参量 把X k 作为输入序列 而输出序列就是x n 把因子WNkn改为WN kn 输入序列的每一个元素除以N 3 4 5FFT IDFT快速算法 IFFT 62 流程如下 共轭法 3 4 5FFT IDFT快速算法 IFFT 63 主题概述 1 绪论2 离散时间系统和离散时间信号的变换3 离散傅里叶变换及其快速计算方法3 1问题的提出3 2DFS 离散傅里叶级数 3 3DFT 有限离散傅里叶变换 3 4FFT 快速离散傅里叶变换 3 5CZT及其快速算法3 6其它变换3 7本章小结4 IIR数字滤波器设计和实现5 FIR数字滤波器设计和实现6 数字信号处理中的有限字长效应 64 DFT不适用于 DFT是均匀分布在Z平面单位圆上N点处的频谱 如果我们取样点不均匀时 则很麻烦 只研究信号的某一频段 要求对该频段取样密集 提高分辨率 如窄带信号的频谱分析 研究非单位圆上的取样值 如频谱峰值探测 需要准确计算N点DFT 且N为大的素数 当x n 是短时间序列时 则得到的频率分辨率2pi N是很低的 提高频谱密度的办法 用补零的方法增加点数 但DFT的点数又大大增加 使计算工作量增大 CZT算法 对z变换采用螺旋线取样 chirp z变换 线性调频z变换 Chirpz transform 3 5线性调频z变换CZT及其快速算法 65 CZT的定义设x n 为已知时间序列 其Z变换的形式为 式中复变量z为 这里s为拉普拉斯变量 是个实数 3 5 1CZT 定义 66 按照上式计算Z x n 必然是从z的实轴开始 为得到任意起始点和以螺旋线规律变化的z值 做如下假设 其中 A W为任意复数 0为A的起始角 第1个取样点 k 0 A0为A起始半径 0为在Z平面中相邻的zk 即zk和zk 1 之间的夹角 W0为任意正数值 M为要分析的点数 不一定定于N 3 5 1CZT 定义 67 CZT表达式 3 5 1CZT 定义 1 取样点沿螺旋线按角度间隔 0分布 0 0时 取样轨迹逆时针旋转 01 随着k的增加 向内盘旋 朝向原点 W0 1 随着k的增加 向外盘旋 3 若W0 1 一段圆弧 若同时A0 1 则为单位圆一部分 此时CZT x n DFT x n M N 68 例3 27A 0 8 exp j pi 6 起始半径0 8 起始角pi 6W 0 985 exp j pi 0 05 W0 1 外旋 取样间隔0 05piM 91 计算点数z1 A W 0 M 1 Zplane z1 3 5 1CZT 定义 69 A 0 8 exp j pi 6 起始半径0 8 起始角pi 6W 1 031 exp j pi 0 05 W0 1 内旋 取样间隔0 05piM 91 计算点数z2 A W 0 M 1 Zplane z2 3 5 1CZT 定义 70 CZT表示为线性卷积形式 利用布鲁斯坦 Bluestein 提出的等式把在CZT x n 中的Wnk因子展开 得 3 5 2CZT 快速算法 计算量 直接计算M点的CZT 需要MN次复数乘法 M N 1 次复数加法 需要快速算法 把上述运算转换为线性卷积形式 从而可以采用FFT算法 提高运算速度 71 3 5 2CZT 快速算法 令 则 式中 72 3 5 2CZT 快速算法 h n x n g n X zk 序列 可以想像为频率随时间 n 线性增长的复指数序列 在雷达中这类信号 称为线性调频信号 chirpSignal 因此 这里z变换称为线性调频z变换 该信号的瞬时频率 i 2 t正比于时间t 因此 该信号的频率随时间线性增长 类比 73 取值范围n 0 N 1 k 0 M 1因果序列g n 的长度为N 而序列h n 的长度为N M 1 并位于区间内 N 1 M 1 所以卷积后的序列y n 的长度为2N M 2 且位于区间内 N 1 N M 2 实际上只要得到区间0 M 1 内的线性卷积结果就能够完成CZT的计算 3 5 2CZT 快速算法 74 用循环卷积代替线性卷积且不产生混迭失真的条件是 循环卷积的点数 周期 应大于或等于2N M 2 但CZT只需要前0 M 1 值 对以后的其它值是否混迭并不关心 因此 可将循环卷积的点数缩减到最小为N M 1 为了基2FFT运算 取L N M 1 L 2m 3 5 2CZT 快速算法 75 g n h n 序列的构造为了进行循环卷积 把h n 以周期L进行周期拓展 再取主值序列 从而得到圆周卷积的一个序列 从n M开始补L N M 1 个零或任意值 补到n L N处 从n L N 1到L 1 取h n 的周期拓展序列h L n 进行循环卷积的另一个序列g n 只需要补零到L即可 3 5 2CZT 快速算法 76 3 5 2CZT 快速算法 h n 的构造 77 CZT快速算法下面说明上述计算方案的具体实现 既然X zk 可用线性卷积表示 我们就可以用DFT来计算线性卷积 只不过DFT换成FFT 从而得到CZT的快速算法 3 5 2CZT 快速算法 L点DFT H k h n 构造 长度M N 1 长度N g n 补零 L点DFT G k G k H k L N M 1 L点IDFT y k 78 3 5 2CZT 快速算法 79 CZT Chirpz transform 变换的Matlab句法为 y czt x m w a y czt x CZT是x信号沿着w和a定义的螺旋线进行的z变换 m是说明了变换的长度 w是z平面上感兴趣的那部分螺旋线上取样点之间的比值 a是螺旋线上的复数起始点 Z平面上的螺旋线 或 线性调频脉冲 定义为 z a w 0 m 1 如果x是矩阵 czt x m w a 是x的列变换 y czt x 使用下列缺省值 m length x w exp j 2 pi m a 1对于这些缺省值 czt返回x信号在单位园上m份等间隔的z变换 也就是x的离散傅立叶变换 或者fft x 空矩阵 说明了这些参数的缺省值 3 5 3CZT Matlab实现 80 例3 28用czt放大一个滤波器频率响应的窄带部分 100 150Hz 首先设计一个滤波器 30thorderLPFfilter cut offfrequencyWn 125 500 use boxcar windowh fir1 30 125 500 boxcar 31 建立频率和CZT初始化参数 Fs 1000 f1 100 f2 150 inHertzm 1024 取样点数w exp j 2 pi f2 f1 m Fs 取样点间隔a exp j 2 pi f1 Fs 取样点起始位置计算滤波器的DFT和CZT变换 y fft h 1000 z czt h m w a 得到频率响应和比较结构 fy 0 length y 1 1000 length y fz 0 length z 1 f2 f1 length z f1 plot fy 1 500 abs y 1 500 axis 150001 2 plot fz abs z axis f1f201 2 3 5 3CZT Matlab实现 81 主题概述 1 绪论2 离散时间系统和离散时间信号的变换3 离散傅里叶变换及其快速计算方法3 1问题的提出3 2DFS 离散傅里叶级数 3 3DFT 有限离散傅里叶变换 3 4FFT 快速离散傅里叶变换 3 5CZT及其快速算法3 6其它变换3 7本章小结4 IIR数字滤波器设计和实现5 FIR数字滤波器设计和实现6 数字信号处理中的有限字长效应 82 二维DFT FFTDFT的局限性和克服方法测不准原理三个局限性短时傅立叶变换STFT时频联合分析滤波器组小波变换框架 Frame 其它离散变换KLT变换Walsh变换Hadamard变换离散余弦变换 3 6其它变换 83 DFT FFT计算不仅仅局限于语音和音乐等一维信号 二维信号的频谱特性也很重要 2DDFT是1DDFT的自然延伸 2DDFT定义如下 2DDFT的可分离性 指首先沿着图像 二维信号 的每一行计算一维变换 然后沿着这一中间结果的每一列再计算一维变换 以此计算二维2D DFT变换 得到最后的结果 3 6 1其它变换 2DDFT FFT 84 对于每个x值 当v 0 1 N 1时 该式为完整的一维傅里叶变换 其中 1DN DFTy v 1DM DFTx u 3 6 1其它变换 2DDFT FFT 85 直接根据前面的公式 计算N点一维傅里叶变换需要N2次复数乘法 N N 1 N2次复数加法运算 快速傅里叶变换只需要约Nlog2N次运算 例如N 1024 直接傅里叶变换需要大约106次操作 快速傅里叶变换只需要104次操作 2D FFT计算量 3 6 1其它变换 2DDFT FFT 86 3 6 2DFT局限性 Heisenberg测不准原理 不定原理 Heisenberg测不准原理 Heisenberg Gabor不定原理 给定信号x t 若 则 当且仅当x t 为高斯信号 即时等号成立 式中 显然 这是方差的标准定义 通常定义2 t 2 分别是信号的时宽和带宽 定义 t 为信号的时宽 带宽积 87 测不准原理是信号处理中的一个重要的基本定理 不可能违背 测不准原理给出了信号时宽 带宽之间的制约关系 对于给定的信号 其时带与带宽的乘积为一个常数 当信号的时宽减少时 其带宽将相应增加 当时宽减到无穷小时 带宽将变成无穷大 如时域的 函数 反之亦然 如频域 处的 函数 时域为ej t cos t jsin t 其在时域的持续时间是 这就是说 信号的时宽与带宽不可能同时趋于无限小 这一基本关系即是我们将要讨论的时间分辨率和频率分辨率的制约关系 在这一基本关系的制约下 人们在竭力探索既能得到好的时间分辨率 或窄的时宽 又能得到好的频率分辨率 或窄的带宽 的信号分析方法 若信号x t 的持续时间是有限的 我们称其为是 紧支撑 CompactSupport 的 其时间的持续区间 如t t1 t2 称为 支撑范围 对频率信号 也使用类似的称呼 3 6 2DFT局限性 Heisenberg测不准原理 88 3 6 2DFT局限性 1807年 傅立叶在论文 OnthePropagationofHeatinSolidBodies 中提出了周期信号傅立叶级数的概念 1822年 傅立叶又在著作 Th orieanalytiquedelachaleurin1822 中提出了非周期信号分解的概念 即傅立叶变换 傅立叶变换不但是重要的数学分支 也是信号分析与信号处理的重要工具 在傅立叶变换应用中 人们早就发现了其不足 Gabor TheoryofCommunication IEEE 1946 体现在三个方面 傅立叶变换在时间和频率定位功能上的局限性 傅立叶变换对于非平稳信号的局限性 傅立叶变换在分辨率上的局限性 89 时间和频率的定位 某一个特定时间t0所对应的频率是多少 或对某一个特定的频率 0所对应的时间是多少 有限长序列的离散付氏变换对定义为 3 6 2DFT局限性 时间和频率的定位 对给定的某一频率k0 为求该频率处的傅氏变换X k0 DFT x n 中的求和范围n需要从0到N 1 即需要x n 整个时域的 知识 反之 如果我们要求出某一时刻n0T处的值x n 由IDFT X k 式 我们仍需要对X k 从0到N 1求和 同样也需要X k 整个频域的 知识 实际上 傅氏变换X k 是信号x n 在整个求和 积分 区间的时间范围内所具有的频率特征的平均表示 反之 也是如此 因此 傅立叶变换不具有时间和频率的 定位 功能 90 例1 设信号x n 由三个不同频率的正弦所组成 xa n 时频分布的三维表示 3 6 2DFT局限性 时间和频率的定位 91 例2 DTMF DualToneMultifrequency 双音多频信号每当按下一个按键时 产生高 低频率的两个音频信号 高 低频音的一个组合表示一个特定的数字0 9 和A D 0123456789 abcd 整个号码的频谱 只有8个频率 3 6 2DFT局限性 时间和频率的定位 92 3 6 2DFT局限性 非平稳信号 非平稳的定义 两类说法 出发点不同 无大碍 非平稳随机信号在平稳和非平稳随机信号定义中 前者是指随机信号的一阶和二阶统计特性 均值 方差 不随时间变换 自相关函数与观察的时间起点无关 若信号是平稳的 则满足维纳 辛钦关系 即功率谱密度与自相关函数互为傅立叶变换 非平稳信号频率随时间变化的信号 称为时变信号 或非平稳信号频率不随时间变化的时不变信号 称为平稳信号注意与随机信号中平稳 非平稳定义的区别 说一个信号是平稳信号 要指明是频率不随时间变化 还是平稳随机信号 若信号是非平稳的 则不满足上述关系 93 例3 一个线性频率调制信号 雷达中的Chirp 为 该信号的瞬时频率 i 2 t正比于时间t 因此该信号的频率是时变的 傅立叶变换对时变信号的局限性 时域波形看 随着时间的增长 信号的振荡 频率 愈来愈快 但从频域波形看不出该信号的频率随时间线性增长的特点 IF曲线也是信号能量的主要集中处 Chirp信号的能量主要集中在时间 频率平面的斜线上 x n 时频分布三维表示 从时频分布能够明确体现出频率与时间的正比关系 3 6 2DFT局限性 时间和频率的定位 94 3 6 2DFT局限性 时间和频率分辨率 分辨率 Resolution 是信号处理的基本概念 它包含了信号的时域和频域两个方面 它是指对信号所能辨别的时域或频域的最小间隔 又称最小分辨细胞 频率分辨率是通过一个频域窗函数观察频谱时所看到的频率的宽度 时间分辨率是通过一个时域窗函数观察信号时所看到的时间的宽度 显然 窗函数越窄 相应的分辨率就越好 分辨能力的好坏取决于 信号的特点 信号的长度 所用的算法 追求时间和频率分辨率同时完好 但由测不准原理可知 时间和频率分辨率不可能同时达到最好 即分辩间隔最小 在实际中 需要根据信号特点和处理需要 选择不同的时间和频率分辨率 对在时域具有瞬变的信号 希望时域的分辨率要好 即时域的观察间隔尽量短 以保证能观察到该瞬变信号发生的时刻及瞬变的形态 当然 这时就忽视频率分辨率 反之 对于时域缓变的信号 就没必要强调时域分辨率 而需好的频率分辨率 好的信号分析算法 应能适应信号的特点而自动调节时间和频率分辨率 95 3 6 2DFT局限性 时间和频率分辨率 理论上 傅立叶变换可以写成如下的内积形式 表示信号和的内积 信号x t 的傅立叶变换等效于x t 和基函数ej t作内积 基函数ej t在频域是位于 处的 函数 因此 当用傅立叶变换来分析信号的频域行为时 它具有最好的频率分辨率 但是 ej t在时域对应的是正弦函数 ej t cos t sin t 其在时域的持续时间是 因此 在时域有着最坏的时间分辨率 对傅立叶反变换 情况正好相反 96 3 6 2DFT局限性 时间和频率分辨率 实际中 将x t 乘以矩形窗进行截短 一个宽度为无穷的矩形窗 即直流信号 的傅立叶变换为 函数 反之亦然 当矩形窗为有限宽时 其傅立叶变换为一Sinc函数 即 矩形窗宽度和其频谱主瓣宽度 T T 成反比 若信号在时域取得越短 即保持高的时间分辨率 那么其频谱的主瓣变宽 必然导致频域的频率分辨率下降 这体现了测不准原理的制约关系 也体现了傅立叶变换中时间和频率分辨率所固有的矛盾 97 DFT的局限性成为寻找新的信号分析和处理的源动力 1932年 Wigner提出了时间 频率联合分布的概念 并将其用于量子物理 1946年 Gabor提出了短时傅立叶变换和Gaber变换概念 从而开始了非平稳信号时频联合分析的研究 20世纪80年代 提出了滤波器组理论 为信号的子带分解提供了有力工具 20世纪80年代后期 发展起来了小波变换 它不仅扩展了信号时频联合分析的概念 而且在信号的分辨率方面具有对信号特点的适应性 为了更一般地讨论信号的分解问题 人们提出了框架 Frame 理论 3 6 2DFT局限性 克服方法 98 3 6 2DFT局限性克服方法 短时傅立叶变换 傅立叶变换中的基函数ej t用下列基函数代替 则短时傅立叶变换STFT或加窗傅立叶变换的定义如下 式中g 是对称的实数窗函数 STFTx t 的意义实际上是用g 沿着t轴滑动 截取一段一段的信号x 然后对其作傅立叶变换 得到 t 平面上二维函数STFTx t g 的作用是保持时域为有限长 一般称作 有限支撑或紧支撑 其宽度越小 则时域分辨率越好 使用不同的基函数可得到不同的分辨率效果 99 可以证明 STFT的基函数gt 在时频平面上具有如下的分辨 细胞 其中心在 t 处 大小为 不管t 取何值 即移到何处 该 细胞 的面积始终保持不变 该面积的大小即是STFT的时频分辨率 3 6 2DFT局限性克服方法 短时傅立叶变换 STFT时间和频率分辨率的自适应性快变信号对应的是高频信号 快变信号需要好的时间分辨率 以观察其快变部分 如尖脉冲等 即观察的时间宽度要小 由测不准原理 时宽 带宽积 知道 该信号频域的分辨率必定要下降 反之 慢变信号对应低频信号 可以降低时域的分辨率 从而在低频处获得好的频率分辨率 希望时频分析算法能自动适应这一要求 显然STFT的 不随t 变化 因而STFT不具备自动调节能力 后面讲的小波变换具备这一能力 100 3 6 2DFT局限性克服方法 时频联合分析 短时傅立叶变换STFT是最简单 最直观的时频联合分析 但STFTx t 中变量t 仍是单独取值 它并不是严格意义上的时频联合分析 20世纪中叶以来 Wigner Ville Gabor和Cohen等陆续给出了一大类真正将时间t和频率 联合起来的分析方法 并在80年代得到广泛应用 Wigner Ville时频分布 Wigner在量子力学中提出了Wigner分布 1948年Ville将它引入信号处理领域 得到了著名的Wigner Ville时频分布 由于x t 在积分中出现了两次 所以又称该式为双线性时频分布 其结果Wx t 是t 的二维函数 它有着一系列好的性质 因此是得到广泛应用的一种信号时频分析方法 101 3 6 2DFT局限性克服方法 时频联合分析 Gabor展开 1946年Gabor提出了信号时频展开的思想 用时间和频率两个坐标同时表示一个信号 即Gabor展开 式中g t 是窗函数 Cm n是展开系数 m代表时域序号 n代表频域序号 Cohen分布 1966年Cohen提出了如下的时频分布 式中g 是处在 平面的权函数 可以证明 若g 1 则Cohen分布即变成Wigner Ville分布 给定不同的权函数 我们可得到不同的时频分布 在80年代前后提出的时频分布有十多种 统称为Cohen类时频分布 简称Cohen类 102 3 6 2DFT局限性克服方法 时频联合分析 对给定的信号x t 人们希望能找到一个二维函数Wx t 它应具有如下基本特性 它是人们最关心的两个物理量时间t和频率 的联合分布函数 它可反映x t 的能量随时间t和频率 变化的形态 既具有好的时间分辨率 同时又具有好的频率分辨率 103 3 6 2DFT局限性克服方法 滤波器组 等带宽的M通道滤波器组 信号的子带分解 信号的多分辨率分析 信号按二进制分解 产生了不均匀的频带分割 从而得到不同的时间和频率分辨率 非均匀树状滤波器组 104 3 6 2DFT局限性克服方法 小波变换定义 在80年代后期及90年代初期所发展起来的小波变换理论已形成了信号分析和信号处理的又一强大的工具 其实 小波分析可看作信号时 频分析的又一种形式 对于给定信号x t 希望找到一个基本函数 t 并通过 t 的伸缩与位移构成一族函数 a b t 式中 a是尺度因子 b是位移 a b均为常数 且a 0 t 又称为基本小波或母小波 a b t 称为小波基函数或小波基 小波变换定义为信号x t 与小波基 a b t 的内积 即 105 3 6 2DFT局限性克服方法 小波变换定义 小波变换对信号的 尺度 位移 联合分析 也是时频分布的一种 尺度 尺度因子a的作用是把母小波 t 作伸缩 由傅立叶变换的性质可知 若 t 的傅立叶变换是 j 则 t a 的傅立叶变换是a ja 若a 1 则 t a 表示将 t 在时间轴上展宽 而将 j 在频率轴上压缩 若a 1 t a 表示将 t 在时间轴上压缩 将 j 在频率轴上展开 若把 t 看成一窗函数 t a 的宽度将随着a的不同而不同 这也同时影响到频域 即 ja 由此可得到不同的时域分辨率和频域分辨率 位移 参数b是沿着时间轴的位移 确定对信号x t 分析的时间位置 也即时间中心 尺度因子a和位移b联合起来确定了对x t 分析的中心位置和分析的时间宽度 106 3 6 2DFT局限性克服方法 小波变换定义 母小波的尺度及位移对分析范围的控制 a 基本小波 b b 0 a 1 c b不变 a 2 d 分析范围 式中的因子是为了保证在不同的尺度a时 小波基 a b t 始终能和母小波 t 有着相同的能量 即 107 3 6 2DFT局限性克服方法 小波变换特点 时间和频率定位功能小波基 a b t 和小波变换WTx a b 在时域中是有限支撑的 具有时间定位功能 a b 和X 在频域中也是有限支撑的 从而具有频率定位功能 小波变换符合时宽 带宽积约束若母小波 t 的时间中心是t0 时宽是 t 则 t a 的时间中心变为at0 时宽变成a t 若母小波 t 的频谱 的频率中心是 0 带宽是 则 t a 的频谱a a 的频率中心变为 0 a 带宽变成 a t a 的时宽 带宽积仍是 t 与尺度a无关 这一方面说明小波变换WTx a b 的时频关系也受测不准原理制约 但另一方面也揭示了小波变换的一个重要性质 即恒Q性质 Q值 品质因数 定义 带宽 中心频率 108 3 6 2DFT局限性克服方法 小波变换特点 小波变换的恒Q性质母小波 t 的Q值 t a 的Q值保持不变 不论a为何值 a 0 t a 始终和 t 具有相同的品质因数Q 恒Q性质是小波变换的一个重要性质 也是区别于其它类型变换且被广泛应用的一个重要原因 下图说明了 和 a 的带宽及中心频率随a变化的情况 109 3 6 2DFT局限性克服方法 小波变换特点 小波变换的时频区间由上页图 当a变小时 对x t 的时域观察范围变窄 但对X 的频率观察范围变宽 且观察的中心频率向高频处移动 反之 当a变大时 对x t 的时域观察范围变宽 频域观察范围变窄 且分析的中心频率向低频处移动 在不同尺度下小波变换所分析的时宽 带宽 时间中心和频率中心的关系 如图 由于恒Q性质 因此在不同尺度下 小波变换的三个时 频分析区间 即三个矩形 的面积保持不变 但提供了一个在时 频平面上长度可调的分析窗口 该分析窗口在高频端的频率分辨率不好 矩形窗的频率边变长 但时域的分辨率变好 矩形的时间边变短 反之 在低频端频率分辨率变好 而时域分辨率变差 但在不同的a值下 分析窗的面积保持不变 也即时频分辨率可以随分析任务的需要作出调整 110 3 6 2DFT局限性克服方法 小波变换特点 小波变换对信号的自适应性小波变换可以自动满足客观实际信号的需要信号中的高频成份往往对应时域中的快变成份 如陡峭的前沿 后沿 尖脉冲等 对这一类信号分析时则要求时域分辨率要好以适应快变成份间隔短的需要 对频域的分辨率则可以放宽 当然 时 频分析窗也应处在高频端的位置 与此相反 低频信号往往是信号中的慢变成份 对这类信号分析时一般希望频率的分辨率要好 而时间的分辨率可以放宽 同时分析的中心频率也应移到低频处 当用较小的a对信号作高频分析时 实际上是用高频小波对信号作时间上的细致观察 当用较大的a对信号作低频分析时 实际上是用低频小波对信号作时间上概貌观察 小波变换的这一特点即既符合对信号作实际分析时的规律 也符合人们的视觉特点 综上所述 由于小波变换具有恒Q性质及自动调节对信号分析的时宽 带宽等一系列突出优点 因此被人们称为信号分析的 数学显微镜 111 3 6 2DFT局限性克服方法 小波变换的容许条件 小波变换的容许条件 小波变换存在的条件 并不是时域的任一函数 t L2 R 都可以充当小波 其可以作为小波的必要条件是其傅里叶变换满足上述容许条件 必有 否则c 必趋于无穷 分母 0 这等效于小波函数 t 是带通函数 由于 因此必有 t 的取值必然是有正有负 也即它是振荡的 上述三条描绘了小波函数的大致特征 即是一带通函数 它的时域波形应是振荡的 此外 从时 频定位的角度 希望是有限支撑的 因此它应是快速衰减的 小波 wavelet 时域有限长 且是振荡的一类函数 112 3 6 2DFT局限性克服方法 小波变换分类 小波变换是80年代后期发展起来的应用数学分支小波理论 法国数学家Meyer Y 地质物理学家Morlet J和理论物理学家Grossman A对小波理论作出了突出的贡献 小波工程应用 法国学者Daubechies I和Mallat S在将小波理论引入工程应用 特别是信号处理领域起到了重要的作用 人们称这些人为 法国学派 小波变换的分类 第一类是所谓的 经典小波 在MATLAB中把它们称作 原始 Crude 小波 这是一批在小波发展历史上比较有名的小波 第二类是Daubecheis构造的正交小波 第三类是由Cohen Daubechies构造的双正交小波 113 3 6 2DFT局限性克服方法 经典小波 Haar小波 源自于数学家Haar于1910年提出的Haar正交函数集 其定义是 Haar小波有很多优点 Haar小波在时域是紧支撑的 Haar小波属正交小波 Haar波是对称的 是目前唯一一个既具有对称性又是有限支撑的正交小波 Haar小波仅取 1和 1 计算简单 Haar小波是不连续小波 使其在实际的信号分析与处理中受到了一定限制 Haar小波有上述优点 在教科书与论文中常被用作范例来讨论 114 3 6 2DFT局限性克服方法 经典小波 Morlet小波 Morlet小波是一个具有高斯包络的单频率复正弦函数 考虑到待分析的信号一般是实信号 所以在MATLAB中将它改造为左式 并取 0 5 该小波不是紧支撑的 t理论上可取 但是当 0 5 或取更大的值时 t 和 在时域和频域都具有很好的集中 Matlab改造 Morl