Matlab数字信号处理

MatlabMatlab数字信号处理数字信号处理11、信号的产生、信号的产生22、信号的运算、信号的运算33、差分方程与、差分方程与ZZ变换变换44、、快速傅里叶变换快速傅里叶变换55、、数字滤波器的设计数字滤波器的设计66、、使用中的一些技巧使用中的一些技巧一、信号的产生一、信号的产生11、、单位采样序列单位采样序列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*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);00.010.020.030.040.050.060.07-1-0.8-0.6-0.4-0.200.20.40.60.81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)。00.010.020.030.040.050.060.07-1-0.8-0.6-0.4-0.200.20.40.60.819、sinc函数信号y=sinc(x);产生周期为2*pi，随x的增加衰减震荡的偶函数，在n*pi处值为零00.010.020.030.040.050.060.07-0.4-0.200.20.40.60.81二、信号的运算二、信号的运算1、信号的延迟给定信号x(n),若信号y1(n)、y2(n)分别定义为：y1(n)=x(n-k)y2(n)=x(n+k)那么，y1(n)是整个x(n)在时间轴上右移k个时间单位所得到的新序列，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、信号的能量及功率信号的能量有如下表示形式：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开始的。例：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]差分方程与差分方程与ZZ变换变换1、离散系统的时域表示由此可见，系统地输出，就是输入与单位抽样响应卷积得到的。例：如下离散系统：系统的单位抽样信号的响应可以通过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=[1zeros(1,63)];%产生单位抽样信号b=[0.20.1];%a=[1–0.4–0.5];h=filter(b,a,x);h1=impz(b,a);figure;stem(h);title(‘Filterfunction’);figure;stem(h1)title(‘Impzfunction’)2115z.04z.011z.02.0)z(h01020304050607000.020.040.060.080.10.120.140.160.180.2Filterfunction05010015000.020.040.060.080.10.120.140.160.180.2Impzfunction2、离散系统的频率响应Matlab中的freqz函数用来计算由a，b构成系统的频率响应。[h,f]=freqz(b,a,n,fs)例：clc;clear;fs=1000;b=[0.20.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’);快速傅里叶变换快速傅里叶变换在matlab中实现fft很简单，只需要通过命令fft来实现，这里需要注意的是采样频率不要太高，否则频谱的信息被掩盖。另外，计算所得的序列中第k点所对应的实际频率为，f=k*fs/N,其中N为进行傅利叶变换的点数，fs为采样频率，一般取信号最大频率的3-5倍。例：clc；clear；loadleleccum;x=leleccum;N=length(x);y=angle(fft(x));fs=100;f=(1:N)*fs/N;plot(f,y);数字滤波器的设计数字滤波器的设计1、IIR数字滤波器的设计MATMB中设汁IIR数字滤波器的步骤总结如下：(1)巴特沃思数字滤波器的设计[b,a]=butter(n,wn);[b,a]=butter(n,wn,’ftype’);……它得到的是一个阶数为n，截止频率为wn的低通滤波器。其中wn是指滤波器的半功率点，取值范围在0-1之间，取1时，为采样频率的一半，如果wn=[w1w2]为两个元素的向量，函数返回的是阶数为2*n的带通滤波器，通带范围为w1-w2,其中ftype代表滤波器的形式，为high时，为高通滤波器，得到阶数为n，截止频率为wn的，高通滤波器。为stop时，得到的是阶数为2*n，阻带为w1-w2的带阻滤波器。例：采样频率为1000Hz的采样信号，设计一个10阶带通butter滤波器，通带范围为100-200Hz，并画出冲击响应曲线。n=5；wn=[100200]/500；[b,a]=butter(n,wn,’bandpass’);[y,t]=impz(b,a,101);stem(t,y);0102030405060708090100-0.25-0.2-0.15-0.1-0.0500.050.10.150.2例:信号采样频率为1000Hz,设计一个阶数为9,截止频率为300Hz的高通巴特沃思滤波器。[b,a]=butter(9,300/500,‘high’);Freqz(b,a,128,1000);滤波器的幅频特性与相频特性如下：050100150200250300350400450500-800-600-400-2000200Frequency(Hz)Phase(degrees)050100150200250300350400450500-400-300-200-1000Frequency(Hz)Magnitude(dB)(2)切比雪夫法设计滤波器切比雪夫法设计滤波器可以分为切比雪夫1法和切比雪夫2法两种，这里我们只介绍切比雪夫1法。其语法结构为：[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=[w1w2]为两个元素的向量，则函数返回的是阶数为2*n的带通滤波器系统函数有理多项式的系数，滤波器的通带范围为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);滤波器的幅频特性与相频特性如下：050100150200250300350400450500-1000-800-600-400-2000Frequency(Hz)Phase(degrees)050100150200250300350400450500-400-300-200-1000Frequency(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,101);stem(t,y);冲激响应曲线如下：0102030405060708090100-0.2-0.15-0.1-0.0500.050.10.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=0.5,B=0.5,C=0布莱克曼窗是三项余弦窗：A=0.42,B=0.5,C=