经典matlab信号处理基础知识.doc

常用函数常用函数 1 1 图形化信号处理工具 图形化信号处理工具 fdatoolfdatool 滤波器设计 滤波器设计 fvtoolfvtool 图形化滤波器参数查看 图形化滤波器参数查看 sptoolsptool 信号处理 信号处理 fvtool b a wintool 窗函数设计 或者使用工具箱 filter design 设计 当使用离散的福利叶变换方法分析频域中的信号时 傅里叶变换时可能引起漏谱 因 此需要采用平滑窗 2 2 数字滤波器和采样频率的关系数字滤波器和采样频率的关系 如果一个数字滤波器的采样率为 FS 那么这个滤波器的分析带宽为 Fs 2 也就是说 这个滤波器只可以分析 0 Fs 2 的信号 举个例字 有两个信号 S1 频率为 20KHz S2 频率为 40KHz 要通过数字方法滤除 S2 你的滤波器的采样率至少要为 Fs 80HKz 否则就分析不到 S2 了 更不可能将它滤 掉了 当然根据采样定理 你的采样率 F0 也必须大于 80HK Fs 和 F0 之间没关系不大 可以任取 只要满足上述关系就行 3 3 两组数据的相关性分析 r corrcoef x y 4 4 expm 求矩阵的整体的 exp 4 4 离散快速傅里叶 fft 信号处理中 傅里叶变换的典型用途是将信号分解成幅值分量和傅里叶变换的典型用途是将信号分解成幅值分量和 频率分量 频率分量 Ft 为连续傅里叶变换 反傅里叶 ifft 5 5 ztrans Z 变换是把离散的数字信号从时域转为频率 6 6 laplace 拉普拉斯变换是把连续的的信号从时域转为频域 7 7 sound x 会在音响里产生 x 所对应的声音 8 8 norm 求范数 det 行列式 rank 求秩 9 9 模拟频率 数字频率 模拟角频率关系模拟频率 数字频率 模拟角频率关系 模拟频率模拟频率 f f 每秒经历多少个周期 单位 每秒经历多少个周期 单位 HzHz 即 即 1 s1 s 模拟角频率模拟角频率 是指每秒经历多少弧度 单位是指每秒经历多少弧度 单位 rad srad s 数字频率数字频率 w w 每个采样点间隔之间的弧度 单位 每个采样点间隔之间的弧度 单位 radrad 2pi f 2pi f w w T T 1010 RMSRMS 求法求法 Rms sqrt sum P 2 或者 norm x sqrt length x varvar 方差的开方是方差的开方是 stdstd 标准差 标准差 RMS 应该是 norm x sqrt length x 吧 求矩阵的 RMSRMS std A std A 1111 ftshiftftshift 作用 将零频点移到频谱的中间作用 将零频点移到频谱的中间 1212 filtfiltfiltfilt零相位滤波零相位滤波 采用两次滤波消除系统的非线性相位 y filtfilt b a x 注意x的长度必须是滤波器阶数的3倍以上 滤波器的 阶数由max length b 1 length a 1 确定 13 h t impz b a n fs 计算滤波器的冲激响应 h为n点冲击响应向量 h x freqz b a n fs 计算频响 有fs时 x为频率f 无fs x为w角频率 常用于查看滤波器的频率特性 14 zplane z p 画图零极点分布图 15 beta unwarp alpha 相位会在穿越 180发生回绕 可将回绕的 1616 stepzstepz 求数字滤波器的阶跃响应求数字滤波器的阶跃响应 h t h t stepz b a n fs stepz b a n fs fvtool b1 a1 b2 a2 bn an fvtool Hd1 Hd2 h fvtool 1515 IIRIIR数字滤波器设计方法数字滤波器设计方法 1 先根据已知带同参数求出最佳滤波器阶数和截止频率 n Wn buttord Wp Ws Rp Rs n Wn buttord Wp Ws Rp Rs s b a butter n Wn ftype s 其中Wp为 0 1之间 Ws为阻带角频率 0 1之间 Rp为通带波纹 或者通带 衰减 Rs为阻带衰减 若果给出的是模拟频率若果给出的是模拟频率fp1fp1通带截止频率通带截止频率 fp2 fp2阻带截止频率阻带截止频率 则则Wp fp1 2 fs Wp fp1 2 fs Ws fp2 2 fs Ws fp2 2 fs 如果给出的是实际数字频率比如如果给出的是实际数字频率比如0 3 pi 0 3 pi 如果给出的是 y filter b a x 或者采用零相位滤波y filtfilt b a x 1515 传统传统 FIRFIR 滤波器滤波器 Ftype 为滤波器类型 比如高通 低通 window 为窗函数类型 Window 窗函数 例子 1 设计一个通带滤波器 带宽为 0 35 0 65 b fir1 48 0 35 0 65 freqz b 1 512 1616窗函数长度窗函数长度 窗函数的长度应等于 窗函数的长度应等于 FIRFIR 滤波器系数个数 滤波器系数个数 即滤波器阶数即滤波器阶数 n 1n 1 1717 加窗函数的加窗函数的 FIRFIR 滤波器长度的确定滤波器长度的确定 17 117 1 buttordbuttord 函数求出最佳滤波器阶数和截止频率 然后用函数求出最佳滤波器阶数和截止频率 然后用 fir1fir1 函数调用 窗函数长函数调用 窗函数长 度为滤波器最度为滤波器最佳阶数佳阶数 n 1n 1 17 2 用窗函数方法设计用窗函数方法设计 FIRFIR 滤波器 由滤波器的过渡带的宽度和选择的窗函数决定滤波器 由滤波器的过渡带的宽度和选择的窗函数决定 这里举一个选用海明窗函数设计低通滤波器的例子 低通滤波器的设计要求是 采样频率为 100Hz 通带截至频率为 3 Hz 阻带截止频率为 5 Hz 通带内最大衰减不高于 0 5 dB 阻带最小衰减不小于 50 dB 使用海明窗函数 确定 N 的步骤有 1 从上表可查得海明窗的精确过渡带宽为 6 6pi N 在有些书中用近似过渡带来计算 这当然没有错 但阶数增大了 相应也增加计算量 2 本低通滤波器的过渡带是 DeltaW Ws Wp 5 3 pi 50 04piDeltaW Ws Wp 5 3 pi 50 04pi 3 N 6 6pi DeltaW 6 6pi 0 04pi 165 所以滤波器 的阶数至少是 165 在该帖子中是用理想低通滤波器的方法来计算的 这里 用 fir1 函数来计算 相应的程序有 fs 100 采样频率 wp 3 pi 50 ws 5 pi 50 deltaw ws wp 过渡带宽 的计算 N ceil 6 6 pi deltaw 1 按海明窗计算所需的滤波器阶数 N0 wdham hamming N 1 海明窗计算 Wn 3 5 100 计算截止频率 b fir1 N Wn wdham H w freqz b 1 db 20 log10 abs H 画频响曲线 plot w fs 2 pi db title 幅度响应 单位 dB grid axis 0 50 100 10 xlabel 频率 单位 Hz ylabel 分贝 set gca XTickMode manual XTick 0 3 5 50 set gca YTickMode manual YTick 50 0 17 数字滤波器函数 Butter cheyshev切比雪夫 b a cheby1 n rp wn options b a besself n wn options b a ellip n rp rs wn options n为阶数 wn为截止频率rad s rs 为阻带起伏 wn在0 1之间 且1对应于采样频率的一半 b a butter n Wn options z p k butter n Wn ftype s z p k butter n Wn ftype A B C D butter n Wn ftype s ftype 对应 high 是高通滤波器的归一化截止频率 low 低通滤波器的归一化截止频率 stop for an order 2 n bandstop digital filter if Wn is a two element vector Wn w1 w2 The stopband is w1 Ws 2 低通滤波器 Wp 和 Ws 为一元矢量且 Wp Ws 3 带通滤波器 Wp 和 Ws 为二元矢量且 WpWs 如 Wp 0 1 0 8 Ws 0 2 0 7 2 2 求出滤波器系数求出滤波器系数 Butter 函数可设计低通 高通 带通和带阻的数字和模拟 IIR 滤波器 其特性为使通带 内的幅度响应最大限度地平坦 但同时损失截止频率处的下降斜度 在期望通带平滑的情 况下 可使用 butter 函数 butter 函数的用法为 b a butter n Wn ftype b a butter n Wn ftype 其中 n 代表滤波器阶数 Wn 代表滤波器的截止频率 二 契比雪夫 I 型 IIR 滤波器的设计 在期望通带下降斜率大的场合 应使用椭圆滤波器或契比雪夫滤波器 在 MATLAB 下可使 用 cheby1 函数设计出契比雪夫 I 型 IIR 滤波器 cheby1 函数可设计低通 高通 带通和带阻契比雪夫 I 型滤 IIR 波器 其通带内为等波 纹 阻带内为单调 契比雪夫 I 型的下降斜度比 II 型大 但其代价是通带内波纹较大 cheby1 函数的用法为 b a cheby1 n Rp Wn ftype 在使用 cheby1 函数设计 IIR 滤波器之前 可使用 cheblord 函数求出滤波器阶数 n 和截止 频率 Wn cheblord 函数可在给定滤波器性能的情况下 选择契比雪夫 I 型滤波器的最小 阶和截止频率 Wn cheblord 函数的用法为 n Wn cheblord Wp Ws Rp Rs 其中 Wp 和 Ws 分别是通带和阻带的拐角频率 截止频率 其取值范围为 0 至 1 之间 当 其值为 1 时代表采样频率的一半 Rp 和 Rs 分别是通带和阻带区的波纹系数 例例 1 1 选择设计 IIR 的 Butterworth 低通滤波器 其 Fs 22050Hz Fp1 3400Hz Fs1 5000Hz Rp 2dB Rs 20dB Fs 22050 Fp1 3400 Fs1 5000 Rp 3 Rs 20 设计指标 wp1 2 Fp1 Fs ws1 2 Fs1 Fs 求归一化频率 确定 butterworth 的最小阶数 N 和频率参数 Wn n Wn buttord wp1 ws1 Rp Rs B A butter N Wn 确定传递函数的分子 分母系数 h f freqz b a Nn Fs value 生成频率响应参数 plot f 20 log abs h 画幅频响应图 plot f angle h 画相频响应图 N Wn buttord Wp Ws Rp Rs 确定 butterworth 的 N 和 Wn N Wn cheblord Wp Ws Rp Rs 确定 Chebyshev 滤波器的 N 和 Wn N Wn cheb2ord Wp Ws Rp Rs 确定 Chebyshev2 滤波器的 N 和 Wn N Wn ellipord Wp Ws Rp Rs 确定椭圆 Ellipse 滤波器 的 N 和 Wn B A butter N Wn type 设计 type 型巴特沃斯 Butterworth 滤波器 filter B A cheby1 N R Wn type 设计 type 型切比雪夫 滤波器 filter B A cheby2 N R Wn type 设计 type 型切比雪夫 滤波器 filter B A ellip N Rp Rs Wn type 设计 type 型椭圆 filter 例子 2 实现了对频率为 20 和 200Hz 单频叠加 cos 信号的低通滤波 使输出仅含有 20Hz 分量 clear fs 1200 采样频率 N 1200 N fs 秒数据 n 0 N 1 t n fs 时间 fL 20 fH 200 sL cos 2 pi fL t sH cos 2 pi fH t s sL sH s in f fft s i 1 250 plot i s in f i figure 1 plot t s title 输入信号 xlabel t s ylabel 幅度 设计低通滤波器 Wp 50 fs Ws 100 fs 截止频率 50Hz 阻带截止频率 100Hz 采样频率 200Hz n Wn buttord Wp Ws 1 50 阻带衰减大于 50db 通带纹波小于 1db 估算得到 Butterworth 低通滤波器的最小阶数 N 和 3dB 截止频率 Wn a b butter n Wn 设计 Butterworth 低通滤波器 h f freqz a b whole fs 求数字低通滤波器的频率响应 f 0 length f 1 fs length f 进行对应的频率转换 figure 2 plot f abs h 绘制 Butterworth 低通滤波器的幅频响应图 title 巴氏低通滤波器 grid sF filter a b s 叠加函数 s 经过低通滤波器以后的新函数 figure 3 plot t sF 绘制叠加函数 S 经过低通滤波器以后的时域图形 xlabel t s ylabel 幅度 SF fft sF N 对叠加函数 S 经过低通滤波器以后的新函数进行 256 点的基 2 快速傅立叶变换 mag abs SF 求幅值 f 0 length SF 1 fs length SF 进行对应的频率转换 figure 4 plot f mag 绘制叠加函数 S 经过低通滤波器以后的频谱图 title 低通滤波后的频谱图 4 4 窗函数法窗函数法 FIRFIR 滤波器设计实验滤波器设计实验 FIR 滤波器的设计任务是选择有限长度的 h n 使传输函数 H 满足技术要求 FIR 滤波 器的设计方法有多种 如窗函数法 频率采样法及其它各种优化设计方法 本实验介绍窗 函数法的 FIR 滤波器设计 窗函数法是使用矩形窗 三角窗 巴特利特窗 汉明窗 汉宁窗和布莱克曼窗等设计出标 准响应的高通 低通 带通和带阻 FIR 滤波器 一 一 firlfirl 函数的使用函数的使用 在 MATLAB 下设计标准响应 FIR 滤波器可使用 firl 函数 firl 函数以经典方法实现加窗 线性相位 FIR 滤波器设计 它可以设计出标准的低通 带通 高通和带阻滤波器 firl 函数的用法为 b firl n Wn ftype Window 各个参数的含义如下 b 滤波器系数 对于一个 n 阶的 FIR 滤波器 其 n 1 个滤波器系数可表示为 b z b 1 b 2 z 1 b n 1 z n n 滤波器阶数 Wn 截止频率 0 Wn 1 Wn 1 对应于采样频率的一半 当设计带通和带阻滤波器时 Wn W1 W2 W1 W2 ftype 当指定 ftype 时 可设计高通和带阻滤波器 Ftype high 时 设计高通 FIR 滤波 器 ftype stop 时设计带阻 FIR 滤波器 低通和带通低通和带通 FIRFIR 滤波器无需输入滤波器无需输入 ftypeftype 参数 参数 Window 窗函数 窗函数的长度应等于 FIR 滤波器系数个数 即阶数 n 1 二 窗函数的使用 在 MATLAB 下 这些窗函数分别为 1 矩形窗 w boxcar n 产生一个 n 点的矩形窗函数 2 三角窗 w triang n 产生一个 n 点的三角窗函数 当 n 为奇数时 三角窗系数为 w k 当 n 为偶数时 三角窗系数为 w k 3 巴特利特窗 w Bartlett n 产生一个 n 点的巴特利特窗函数 巴特利特窗系数为 w k 巴特利特窗与三角窗非常相似 巴特利特窗在取样点 1 和 n 上总以零结束 而三角窗在这 些点上并不为零 实际上 当 n 为奇数时 bartlett n 的中心 n 2 个点等效于 triang n 2 4 汉明窗 w hamming n 产生一个 n 点的汉明窗函数 汉明窗系数为 w k 1 0 54 0 46cos k 0 n 1 5 汉宁窗 w hanning n 产生一个 n 点的汉宁窗函数 汉宁窗系数为 w k 0 5 1 cos k 1 n 6 布莱克曼窗 w Blackman n 产生一个 n 点的布莱克曼窗函数 布莱克曼窗系数为 w k 0 42 0 5cos 2 0 8cos 4 k 1 n 与等长度的汉明窗和汉宁窗相比 布莱克曼窗的主瓣稍宽 旁瓣稍低 7 凯泽窗 w Kaiser n beta 产生一个 n 点的凯泽窗数 其中 beta 为影响窗函数旁瓣 的 参数 其最小的旁瓣抑制 与 的关系为 0 1102 0 87 50 0 5842 21 0 4 0 07886 21 21 50 0 0 8 output1 i 0 elseif corr 0 1 output1 i 1 end end for i 1 count4 n 160 i 1 1 160 160 i 1 if output1 i 1 switch abs recoefs1 i case abs recoefs1 i 0 002 elseif corr 0 1 output2 i 1 end end for i 1 count5 n 160 i 1 1 160 160 i 1 if output2 i 1 switch abs recoefs2 i case abs recoefs2 i 0 002 x noissin figure subplot 611 plot x ylabel x 使用 db5 对信号进行 5 层分解 c l wavedec x 5 db5 for i 1 5 对分解的第 5 层到第 1 层的低频系数进行重构 a wrcoef a c l db5 6 i subplot 6 1 i 1 plot a ylabel a num2str 6 i end figure subplot 611 plot x ylabel x for i 1 5 对分解的第 5 层到第 1 层的高频系数进行重构 d wrcoef d c l db5 6 i subplot 6 1 i 1 plot d ylabel d num2str 6 i end 小波分解的细节信号是由白噪声分解得到的 而正弦信号可以在图中的近似信号 a5 得到 因为这一层噪声对正弦信号的影响可以忽略了 2 正弦信号加三角波正弦信号加三角波 下面通过使用小波分析一个由正弦信号 正弦信号的周期为 20 加三角波组成的信号 说明小波分析如何分离这两种信号 程序代码如下 load trsin x trsin subplot 711 plot x axis tight ylabel x 使用 db5 对信号进行 6 层分解 c l wavedec x 6 db5 for i 1 6 对分解的第 6 层到第 1 层的低频系数进行重构 a wrcoef a c l db5 7 i subplot 7 1 i 1 plot a axis tight ylabel a num2str 7 i end figure subplot 711 plot x axis tight ylabel x for i 1 6 对分解的第 5 层到第 1 层的高频系数进行重构 d wrcoef d c l db5 7 i subplot 7 1 i 1 plot d axis tight ylabel d num2str 7 i end 例子例子 3 应用小波分析来识别某一频率区间上的信号应用小波分析来识别某一频率区间上的信号 在本例中 使用小波分析一个由三个不同频率的正弦信号叠加的信号 看是否能将这三 个正弦信号区分开来 结果证明小波分析可以很好地识别某一频率区间的信号 程序代码如下 load sumsin x sumsin figure subplot 611 plot x ylabel x title 原始信号以及各层近似信号 使用 db3 小波进行 5 层分解 c l wavedec x 5 db3 for i 1 5 对分解的第 5 层到第 1 层的低频系数分别进行重构 a wrcoef a c l db3 6 i subplot 6 1 i 1 plot a ylabel a num2str 6 i end figure subplot 611 plot x ylabel x for i 1 5 对分解的第 5 层到第 1 层的高频系数进行重构 d wrcoef d c l db3 6 i subplot 6 1 i 1 plot d ylabel d num2str 6 i end 在本例中 该信号是由周期分别为 200 20 2 的信号组成的 它们的采样周期均 为 1 为方便起见 在此分别称为低频 中频和高频的正弦信号 从图中可以看出 低 频 中频和高频信号分别对应于分解的近似信号 a4 细节信号 d4 以及细节信号 d1 matlabmatlab 小波通用函数小波通用函数 Allnodes 计算树结点 appcoef 提取一维小波变换低频系数 appcoef2 提取二维小波分解低频系数 bestlevt 计算完整最佳小波包树 besttree 计算最佳 优 树 biorfilt 双正交样条小波滤波器组 biorwavf 双正交样条小波滤波器 centfrq 求小波中心频率 cgauwavf Complex Gaussian 小波 cmorwavf coiflets 小波滤波器 cwt 一维连续小波变换 dbaux Daubechies 小波滤波器计算 dbwavf Daubechies 小波滤波器 dbwavf W W dbN N 1 2 3 50 ddencmp 获取默认值阈值 软或硬 熵标准 depo2ind 将深度 位置结点形式转化成索引结点形式 detcoef 提取一维小波变换高频系数 detcoef2 提取二维小波分解高频系数 disp 显示文本或矩阵 drawtree 画小波包分解树 GUI dtree 构造 DTREE 类 dwt 单尺度一维离散小波变换 dwt2 单尺度二维离散小波变换 dwtmode 离散小波变换拓展模式 dyaddown 二元取样 dyadup 二元插值 entrupd 更新小波包的熵值 fbspwavf B 样条小波 gauswavf Gaussian 小波 get 获取对象属性值 idwt 单尺度一维离散小波逆变换 idwt2 单尺度二维离散小波逆变换 ind2depo 将索引结点形式转化成深度 位置结点形式 intwave 积分小波数 isnode 判断结点是否存在 istnode 判断结点是否是终结点并返回排列值 iswt 一维逆 SWT Stationary Wavelet Transform 变换 iswt2 二维逆 SWT 变换 leaves Determine terminal nodes mexihat 墨西哥帽小波 meyer Meyer 小波 meyeraux Meyer 小波辅助函数 morlet Morlet 小波 nodease 计算上溯结点 nodedesc 计算下溯结点 子结点 nodejoin 重组结点 nodepar 寻找父结点 nodesplt 分割 分解 结点 noleaves Determine nonterminal nodes ntnode Number of terminal nodes ntree Constructor for the class NTREE orthfilt 正交小波滤波器组 plot 绘制向量或矩阵的图形 qmf 镜像二次滤波器 rbiowavf Reverse biorthogonal spline wavelet filters read 读取二进制数据 readtree 读取小波包分解树 scal2frq Scale to frequency set shanwavf Shannon wavelets swt 一维 SWT Stationary Wavelet Transform 变换 swt2 二维 SWT 变换 symaux Symlet wavelet filter computation symwavf Symlets 小波滤波器 thselect 信号消噪的阈值选择 thodes References treedpth 求树的深度 treeord 求树结构的叉数 upcoef 一维小波分解系数的直接重构 upcoef2 二维小