1、Matlab Matlab 数字信号处理 数字信号处理 1 1 、信号的产生 、信号的产生 2 2 、信号的运算 、信号的运算 3 3 、差分方程与 、差分方程与 Z Z 变换 变换 4 4 、 、 快速傅里叶变换 快速傅里叶变换 5 5 、 、 数字滤波器的设计 数字滤波器的设计 6 6 、 、 使用中的一些技巧 使用中的一些技巧一、信号的产生 一、信号的产生 1 1 、 、 单位采样序列 单位采样序列 x=zeros(1,n); x=zeros(1,n); x(k)=1; x(k)=1;2 、单位阶跃序列 x=ones(1,n); 3 、正弦序列 n=0:N-1; x=sin(2*pi*f
2、*n*Ts+fai); 、复正弦序列 n=0:N-1; x=exp(j*w*n);5 、指数序列 n=1:N; x=a.n; % 此处必须用. 而不能直接用 6 、随机序列 rand(m,n) % 产生m行，n 列的在0 ，1 上服 %从均匀分布的随机数矩阵 randn(m,n) % 产生均值为0 ，方差为1 的高斯随机 % 序列 7 、方波信号 x=square(t,duty); 产生周期为2*pi ，幅值为正负1 的方波信号，其中 duty为正幅值部分占周期的百分数例： t = 0:.0001:.0625; y = square(2*pi*30*t); plot(t,y); 0 0.01
3、0.02 0.03 0.04 0.05 0.06 0.07 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 8 、三角波( 锯齿波) sawtooth(t,width); 产生周期为2*pi 幅值为正负1 的三角波， width为宽度，取0-1 之间的数 例： t = 0:.0001:.0625; y = sawtooth(2*pi*30*t,1); plot(t,y) ； sawtooth函数类似于sin函数，其中width用于调 整三角波峰值位置，sawtooth(t,1) 等价于 sawtooth(t) 。0 0.01 0.02 0.03 0.04
4、0.05 0.06 0.07 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 9 、sinc 函数信号 y=sinc(x); 产生周期为2*pi ，随x 的增加衰减震荡的 偶函数，在n*pi 处值为零 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1二、信号的运算 二、信号的运算 1 、信号的延迟 给定信号x(n), 若信号y1(n)、y2(n)分别定义为： y1(n)=x(n-k) y2(n)=x(n+k) 那么，y1(n)是整个x(n) 在时间轴上右移k 个时间 单位所
5、得到的新序列， y2(n) 是整个x(n) 在时间 轴上左移k 个时间单位所得到的结果。 编程实现： function y,n=sig_shift(x,m,n0) m 为输入x 的下标；n0 为延迟单位 n=m+n0 ； y=x ； 2 、相加、相乘 x(n)=x1(n)+x2(n); x(n)=x1(n)*x2(n) 当两个向量相乘时，若用.* 表示数组相乘， 此时，x1中对应元素与x2 中对应元素相乘，所 得结果作为结果数组( 矩阵) ，要求两原始数组 中元素个数相同，如果采用* 是进行向量( 矩阵) 的乘法，相加时要求两原始数组中元素个数相 同。3 、信号的能量及功率 信号的能量有如下表
6、示形式： E=sum(abs(x).2); 信号的功率有如下表示形式： E=sum(abs(x).2)/length(x); 4 、信号的折叠 信号折叠就是对x(n)每一项对n=0 的纵坐标 进行折叠, 即: y(n)=x(-n) y(n) 与x(n)关于n=0 对称; y=fliplr(x); n=-fliplr(n); 在实际应用中， fliplr 的主要作用是把序列倒 转， 例：x=1,2,3;4,5,6; y=fliplr(x); %y=3,2,1;6,5,4 6 、信号的卷积 Matlab提供了内部函数conv来实现两个有限长 序列的卷积，该函数假定两个序列的是从n=0 开始的。 例
7、：x=3,11,7,0,-1,4,2; h=2,3,0,-5,2,1; y=conv(x,h); %y=6,31,47,6,-51,-5,41,18,-22,-3,8,2 ( 共n+m-1 项) 6 、信号的相关 (1) 两个序列x(n) 和y(n) 的相关可以看作是 x(n) 与y(-n) 的卷积。同理，信号x(n) 的自 相关即为x(n) 与x(-n) 的卷积。 (2)xcorr(x) 或xcorr(x,y) 例:t=1:5; y=xcorr(t); %y=5,14,6,40,55,40,26,14,5差分方程与 差分方程与 Z Z 变换 变换 1 、离散系统的时域表示 由此可见，系统地输
8、出，就是输入与单位抽样响应卷 积得到的。 例：如下离散系统： 系统的单位抽样信号的响应可以通过filter 函数和impz 函数实现。 (1)filter 函数 因为一个离散系统 可以看作是一个滤波器，该函数 就是利用滤波器来实现的。 它有如下两种形式： y=filter(b,a,x) 由上图可以知道：b=0.2 ，0.1 a=1 ，-0.4 ，-0.5 则系统的单位抽样响应为：h=filter(b,a,x) (2)impz 函数实现 Impz(b,a) 可以直接得到单位抽样响应，并画出 响应图形 程序实现： clc; clear; x=1 zeros(1,63); % 产生单位抽样信号 b=
9、0.2 0.1; % a=1 0.4 0.5; h=filter(b,a,x); h1=impz(b,a); figure; stem(h); title(Filter function); figure; stem(h1) title(Impz function) 2 1 1 5z . 0 4z . 0 1 1z . 0 2 . 0 ) z ( h 0 10 20 30 40 50 60 70 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Filter function 0 50 100 150 0 0.02 0.04 0.06 0.08
10、 0.1 0.12 0.14 0.16 0.18 0.2 Impz function 2 、离散系统的频率响应 Matlab中的freqz 函数用来计算由a ，b 构成系统 的频率响应。 h,f=freqz(b,a,n,fs) 例： clc; clear; fs=1000; b=0.2 0.1; a=1 0.4 0.5; h,f=freqz(b,a,256,fs); mag=abs(h); ph=angle(h); ph=ph*180/pi; figure; plot(f,mag); title(f-m); figure; plot(f,ph); title(f-p);快速傅里叶变换 快速傅里
11、叶变换 在matlab 中实现fft 很简单，只需要通过命 令fft 来实现，这里需要注意的是采样频 率不要太高，否则频谱的信息被掩盖。 另外，计算所得的序列中第k 点所对应的 实际频率为，f=k*fs/N, 其中N 为进行傅利 叶变换的点数，fs 为采样频率，一般取信 号最大频率的3-5 倍。 例： clc ； clear ； load leleccum; x=leleccum; N=length(x); y=angle(fft(x); fs=100; f= (1:N)*fs/N; plot(f,y);数字滤波器的设计 数字滤波器的设计 1 、IIR 数字滤波器的设计 MATMB中设汁IIR
12、 数字滤波器的步骤总 结如下： (1)巴特沃思数字滤波器的设计 b,a=butter(n,wn); b,a=butter(n,wn,ftype); 它得到的是一个阶数为n ，截止频率为wn 的低通滤波器。其中wn 是指滤波器的半功率 点，取值范围在0-1 之间，取1 时，为采样频率 的一半，如果wn=w1 w2 为两个元素的向量， 函数返回的是阶数为2*n 的带通滤波器，通带 范围为w1-w2, 其中ftype 代表滤波器的形式，为 high时，为高通滤波器，得到阶数为n ，截止频 率为wn 的，高通滤波器。为stop 时，得到的是 阶数为2*n ，阻带为w1-w2 的带阻滤波器。 例：采样频
13、率为1000Hz的采样信号，设 计一个10阶带通butter滤波器，通带范 围为100-200Hz，并画出冲击响应曲线。 n=5 ； wn=100 200/500 ； b,a=butter(n,wn,bandpass); y,t=impz(b,a,101); stem(t,y);0 10 20 30 40 50 60 70 80 90 100 -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 例: 信号采样频率为1000Hz, 设计一个阶数为 9, 截止频率为300Hz 的高通巴特沃思滤波 器。 b,a=butter(9,300/500,high)
14、; Freqz(b,a,128,1000); 滤波器的幅频特性与相频特性如下：0 50 100 150 200 250 300 350 400 450 500 -800 -600 -400 -200 0 200 Frequency (Hz) Phase (degrees) 0 50 100 150 200 250 300 350 400 450 500 -400 -300 -200 -100 0 Frequency (Hz) Magnitude (dB) (2)切比雪夫法设计滤波器 切比雪夫法设计滤波器可以分为切比雪夫1 法 和切比雪夫2 法两种，这里我们只介绍切比雪 夫1 法。 其语法结构为
15、： b,a=cheby1(n,Rp,wn); b,a=cheby1(n,Rp,wn,ftype); 设计的是一个阶数为n ，截止频率为wn ，通带 波纹衰减为Rp 的低通滤波器。返回值a,b 分别是 阶数为n+1 的向量，表示滤波器系统函数的分 母和分子的多项式系数，滤波器的传递函数可 以表示为： 函数的截止频率wn 是指通带的边缘，在那滤波 器的幅度响应为-RpdB ，wn 的取值范围为0-1 ， 其中1 表示采样频率的一半。越小的通带波 纹，会导致越大的过渡带宽。如果wn=w1 w2 为两个元素的向量，则函数返回的是阶数为 2*n 的带通滤波器系统函数有理多项式的系 数，滤波器的通带范围为
16、w1-w2 。 用b,a=cheby1(n,Rp,wn,ftype) 设计搞通和带阻 滤波器，ftype 确定滤波器的形式，为high时， 为阶数为n ，截至频率为wn 的高通滤波器；为 stop时，阶数为2*n ，阻带为w1-w2 的带阻滤波 器。 例: 信号采样频率为1000Hz, 设计一个阶数为 9, 截止频率为300Hz 的低通切比雪夫滤波 器，其中滤波器在通带的波纹为0.5dB 。 b,a=cheby1(9,0.5,300/500); Freqz(b,a,512,1000); 滤波器的幅频特性与相频特性如下：0 50 100 150 200 250 300 350 400 450 5
17、00 -1000 -800 -600 -400 -200 0 Frequency (Hz) Phase (degrees) 0 50 100 150 200 250 300 350 400 450 500 -400 -300 -200 -100 0 Frequency (Hz) Magnitude (dB) 例：采样频率为1000Hz的采样信号，设计一个 10阶带通切比雪夫滤波器，通带范围为100- 200Hz，滤波器在通带的波纹为0.5dB ，并画出 冲激响应曲线。 n=10; Rp=0.5; wn=100,200/500; b,a=cheby1(n,Rp,wn); y,t=impz(b,a
18、,101); stem(t,y); 冲激响应曲线如下：0 10 20 30 40 50 60 70 80 90 100 -0.2 -0.15 -0.1 -0.05 0 0.05 0.1 0.15窗函数 窗函数布莱克曼窗，海明窗，汉宁窗，以及矩形窗都是广义 余弦窗的特殊情形。这些窗可以看作是频率为0， 2pi/(N-1) 和4pi/(N-1) 的余弦序列的线性组合。N 代表创 的长度。此类窗的生成方法如下： Ind=(0:n-1)*2*pi/(n-1) W=A-B*cos(ind)+C*cos(2*ind) 海明窗和汉宁窗是两项余弦窗，对海明窗： A=0.54,B=0.46,C=0 对汉宁窗：A
19、=0.5,B=0.5,C=0 布莱克曼窗是三项余弦窗：A=0.42,B=0.5,C=0.08 三种窗可通过以下三个函数实现： hamming ，hanning ，blackman使用中的一些技巧 使用中的一些技巧 1 、数据的装入与存储 save ；load ； 2 、子函数的编写 3 、数据的保存与读取 4 、路径设置 5 、已有程序参考界面制作 界面制作 参考书目：精通matlab综合辅导与指南 菜单 控制框 函数回调的考虑 指针和鼠标按钮事件 中断回调的规则 举例 对话框和请求程序菜单 菜单 主要包括建立菜单，建立子菜单，去除默认菜单，设 置菜单回调函数 clear; clc; hf=f
20、igure; % 建立窗口并返回句柄hf set(hf,menubar,none) % 去除系统菜单 hm=uimenu(hf,label, 我的菜单); hm_exgrid=uimenu(hm,label, 调用内联函数 ,callback,alnn) % 建立子菜单 调用回调函数 ，并调用函数alnn（自己 定义） function alnn() t=linspace(0,1,1024); f=1; x=sin(2*pi*f*t); y=cos(2*pi*f*t) z=x*y; surf(t,t,z);录制语音信号的程序 录制语音信号的程序 fs=11025; % 采样频率 duratio
21、n=2; % 录音时间 fprintf(按任意键开始 %g 秒录音：, duration); pause fprintf(录音中.); y=wavrecord(duration*fs, fs); % duration*fs 是录音资料点数 fprintf(录音结束n); fprintf(按任意键后开始播放：); pause wavplay(y,fs);设计实例 设计实例 1 1 1 、建立一个*.m 文件，实现信号发生功 能(50Hz 正弦) 2 、对该信号进行FFT 3 、将上述程序转化到子函数中 4 、利用按钮调用这两个功能函数，并绘 制图形输出 %signal1.m %= %产生一个频率
22、为50Hz 的正弦信号，取样时间 t=00.2s, 取样频率500Hz clear; clc; t=linspace(0,0.2,100); f=50 %fs=320 x=sin(2*pi*f*t); plot(t,x); sinwave=x; save sinwave; %= %fftrans.m %= %提取刚才产生的sinwave 信号，对其进行fft ， 已知信号采样频率500Hz clear; clc; load sinwave x=sinwave; %fs=320 N=length(x); y=abs(fft(x); f=(1:N)*fs/N plot(f,y); %usefunc
23、tion.m %= %编制子函数，产生制定频率，幅度，相位的 正弦信号，取样时间t=00.2s, 取样频率为信号 频率的10 倍 function y=usefunction1(A,f,fai) fs=10*f; t=linspace(0,0.2,0.2*fs); x=A*sin(2*pi*f*t+fai); plot(t,x); sinwave_usefunction=x; save sinwave_usefunction; save fs; %fft_usefunction.m %= function y=usefunction2(x,fs) N=length(x); y=abs(fft(
24、x); f=(1:N)*fs/N plot(f,y); fft_usefunction=y; save fft_usefunction; %sin_menu.m %= % 建立一个菜单，其有两个子菜单，分别调用上述两个程序 clear; clc; hf=figure; % 建立窗口并返回句柄hf set(hf,menubar,none) % 去除系统菜单 hm=uimenu(hf,label, 我的菜单); hm_exgrid=uimenu(hm,label, 正弦信 号,callback,usefunction1(2,50,0) load sinwave_usefunction x= sinwave_usefunction load fs hm_exgrid=uimenu(hm,label,fft,callback, y=usefunction2(x,fs) )
