MATLAB程序设计 课件 第11章 信号处理

1MATLAB程序设计《MATLAB程序设计》2第十一章信号处理11.1MATLAB信号处理基础知识 11.2统计信号处理 11.3IIR滤波器 11.4FIR滤波器 11.5特殊波形发生函数 11.1MATLAB信号处理基础知识《MATLAB程序设计》3信号采样即为对连续信号离散化过程，假设以频率fs对信号进行采样，则其采样间隔Δt可表示为11.1.1信号采样

假设对一个连续的正弦信号进行离散采样，令正弦信号表达式为采用fs对上述信号进行采样，则对应某个时刻t，有得到离散化后的信号为11.1MATLAB信号处理基础知识《MATLAB程序设计》4例11-1假设一正弦信号幅值A为1，频率f为50Hz，相位φ为20°，常偏量C为0，t取值范围为0~0.1。令采样频率fs为500Hz，请对该信号进行离散化，并绘制采样后的离散信号序列。11.1.1信号采样

A=1;f=50;fs=500;detaT=1/fs;t1=0;t2=0.1;t=t1:detaT:t2;phi=20/180*pi;%信号采样模拟y=A*sin(2*pi*f*t+phi);%信号采样点显示figureplot(t，y，'k--o'，'MarkerFaceColor'，'k');xlabel('t’);

ylabel('y');11.1MATLAB信号处理基础知识《MATLAB程序设计》5工程实际中采样获得的信号，均会受到不同程度的环境噪声和采集系统噪声影响。通常将信号外部噪声干扰看作高斯白噪声。11.1.2噪声模拟

MATLAB函数包normrnd可以产生指定均值和标准差的高斯白噪声，其调用格式如下：r=normrnd(mu，sigma)r=normrnd(mu，sigma，sz1，…，szN)r=normrnd(mu，sigma，sz)11.1MATLAB信号处理基础知识《MATLAB程序设计》6例11-2令一信号如例11-1所示，采样频率fs为1500Hz，在采集过程中，考虑环境和采集系统噪声影响，令噪声符合均值为0，标准差为0.5的高斯白噪声分布，请仿真模拟受噪声干扰的采样信号。11.1.2噪声模拟

A=1;f=50;fs=1500;detaT=1/fs;t1=0;t2=0.1;t=t1:detaT:t2;phi=20/180*pi;%信号采样模拟y=A*sin(2*pi*f*t+phi);%噪声模拟noise=normrnd(0，0.2，size(y));%噪声干扰后的信号yNoise=y+noise;%信号采样点显示figureplot(t，yNoise，'r-'，'MarkerFaceColor'，'k');holdonplot(t，y，'k--')xlabel('t');ylabel('y');legend('yNosie'，'y')11.1MATLAB信号处理基础知识《MATLAB程序设计》7信号频谱分析是进行信号处理的常规手段，通过傅里叶变换，将时域信号转变为频域信号，从频谱分布的角度对信号进行分析，有利于获得信号隐藏的更为丰富的信息，并进一步用于指导工程实践。11.1.3频谱分析

假设一个离散信号序列为xi（i=0，1，2，…，N-1）则对xi进行傅里叶变换为

其中exp(·)表示以e为底的指数运算，k表示数值频率，N表示离散数据点数，Xk表示对应数值k的傅里叶变换结果。11.1MATLAB信号处理基础知识《MATLAB程序设计》8傅里叶逆变换为11.1.3频谱分析

通过傅里叶逆变换可将频域信息再次转换为时域信息。上述傅里叶变换和逆变换是数学上的定义。在工程实际中，为了赋予数值频率k以物理意义，通常将数值频率k转换为物理频率f(单位：Hz)，转换过程如下，11.1MATLAB信号处理基础知识《MATLAB程序设计》911.1.3频谱分析

在MATLAB中，通常采用快速傅里叶算法实现信号时域和频域之间的转换，对应的函数包为fft，其具体调用格式如下，Y=fft(x)Y=fft(x，n)Y=fft(x，n，dim)例11-3假设信号符合下式其中，A1为0.7，A2为1，f1和f2分别为50Hz和120Hz，采用1kHz采样频率对上述信号进行采样，采样时长为1.5s，noise符合均值为0，标准差为2的高斯白噪声随机分布，使用傅里叶变换求噪声中隐藏的信号频率分量11.1MATLAB信号处理基础知识《MATLAB程序设计》1011.1.3频谱分析

fs=1e3;detaT=1/fs;t1=0;t2=1.5;t=t1:detaT:t2;A1=0.7;A2=1;f1=50;f2=120;x=A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+normrnd(0，2，size(t));%信号显示figureplot(t，x);xlabel('t');ylabel('y');%傅里叶变换L=length(x);Xk=fft(x);P2=abs(Xk/L);P1=P2(1:L/2+1);P1(2:end-1)=2*P1(2:end-1);f=fs*(0:(L/2))/L;%频谱显示figureplot(f，P1);xlabel('f(Hz)')ylabel('|Xk|')11.2统计信号处理《MATLAB程序设计》1111.2.1自相关计算

离散序列自相关定义如下

其中m表示序列平移量，n表示数据序列号，r表示相关性，上式描述了序列x(n)与平移m后的x(n-m)之间的相似性。通常情况下会将相关量r进行归一化，其表达式为其中ρ称为相关系数11.2统计信号处理《MATLAB程序设计》1211.2.1自相关计算

MATLAB中的xcorr()函数可实现离散序列自相关计算。xcorr()函数调用格式如下：r=xcorr(x，y)r=xcorr(x)r=xcorr(_，maxlag)r=xcorr(_，scaleopt)[r，lags]=xcorr(_)11.2统计信号处理《MATLAB程序设计》1311.2.1自相关计算

例11-4假设离散序列取值表达式为，n取值分别为0，1，2，…，15，计算x离散序列自相关系数，并绘制自相关系数ρ与平移量m之间的关系。n=0:15;x=0.84.^n;[r，m]=xcorr(x，"normalized");figurestem(m，r，'k');xlabel('m’)

ylabel('r')11.2统计信号处理《MATLAB程序设计》1411.2.2互相关计算

离散序列互相关定义如下此时r描述了序列x(n)与序列y(n)平移m后的相似性。互相关量r进行归一化，其表达式为其中ρ称为互相关系数。11.2统计信号处理《MATLAB程序设计》1511.2.2互相关计算

例11-5假设x和y序列取值满足下式其中A取值为1，fx和fy分别取值为50和100，t取值范围为0~0.02，取值间隔为0.001，试计算x与y之间互相关系数ρ，并绘制ρ与偏移量m之间的关系图11.2统计信号处理《MATLAB程序设计》1611.2.2互相关计算

detaT=0.001;t0=0;t1=0.02;t=t0:detaT:t1;A=1;fx=50;fy=100;x=A*sin(2*pi*fx*t);y=A*cos(2*pi*fy*t);%互相关计算[r，m]=xcorr(x，y，"normalized");%结果显示figurestem(m，r，'k');xlabel('m')ylabel('r')11.2统计信号处理《MATLAB程序设计》1711.2.3功率谱

功率谱估计被广泛应用于分析平稳各态遍历随机信号频率成分，并且已被成功应用到雷达信号处理、机械故障诊断等工程实际中。给定一个标准的正弦信号，可以通过傅里叶变换来分析它的频率成分。然而，在实际工程应用中，由于存在各种噪声干扰，我们得到的信号往往不是理想的。对于这种存在噪声干扰的信号，可以通过功率谱对信号的频率成分进行分析，得到一个“非精确”的功率谱来对真实随机信号的功率谱进行估计。11.2统计信号处理《MATLAB程序设计》1811.2.3功率谱

例11-6假设随机信号满足下式

令A1和A2均为1，f1为50Hz，f2为135Hz，φ1和φ2分别为0°和45°，Noise符合均值为0，标准差为0.2的高斯白噪声分布。假设采样频率fs为500Hz，请采用周期图法对该随机信号进行功率谱分析。11.2统计信号处理《MATLAB程序设计》1911.2.3功率谱

A1=10;A2=10;f1=50;f2=135;fs=500;T=1/fs;t1=0;t2=1;t=t1:T:t2;phi1=0;phi2=45/180*pi;y=A1*sin(2*pi*f1*t+phi1)+A2*sin(2*pi*f2*t+phi2)+normrnd(0，0.2，1，length(t));%信号波形显示figureplot(t，y);xlim([0，0.2]);%显示一部分波形xlabel('t');ylabel('y');%功率谱分析figureperiodogram(y，[]，length(y)，fs);11.3IIR滤波器《MATLAB程序设计》2011.3.1IIR滤波基本原理

IIR滤波器是一种数字滤波器，其系统函数为

其中，H(z)表示系统函数，X(z)表示系统输入函数，Y(z)表示系统输出函数，z为时延符号，ak和bk表示滤波系数。11.3IIR滤波器《MATLAB程序设计》2111.3.2巴特沃斯滤波器

butter函数用于设计butterworth数字滤波器，其调用格式如下：[b，a]=butter(n，Wn)[b，a]=butter(n，Wn，ftype)[z，p，k]=butter(___)[A，B，C，D]=butter(___)[___]=butter(___，'s')11.3IIR滤波器《MATLAB程序设计》2211.3.2巴特沃斯滤波器

例11-7假设采样频率为2000Hz，设计一个8阶高通Butterworth滤波器，其中截止频率为300Hz。fs=2000;n=8;Wn1=300;Wn=Wn1/(fs/2);[b,a]=butter(n,Wn,'high');%频率响应sampleN=128;freqz(b,a,sampleN,fs);axis([0500-400100]);11.3IIR滤波器《MATLAB程序设计》2311.3.2巴特沃斯滤波器

例11-8假设处理的数字信号采样频率为2000Hz，设计一个10阶的带通Butterworth数字滤波器，通带为150Hz到300Hz，并绘制其脉冲响应曲线。fs=2000;n=10;W1=150;W2=300;sampleN=101;Wn=[W1,W2]/(fs/2);[b,a]=butter(n,Wn);[y,t]=impz(b,a,sampleN);stem(t,y,'k');xlabel('t');ylabel('y');11.3IIR滤波器《MATLAB程序设计》2411.3.3切比雪夫滤波器

cheby1函数cheby1函数用于设计切比雪夫数字滤波器中的chebyshevI型。其调用格式为：[b，a]=cheby1(n，Rp，Wp)[b，a]=cheby1(n，Rp，Wp，ftype)[z，p，k]=cheby1(___)[A，B，C，D]=cheby1(___)[___]=cheby1(___，'s')本节将重点介绍使用cheby1和cheby2设计切比雪夫滤波器的方法。11.3IIR滤波器《MATLAB程序设计》2511.3.3切比雪夫滤波器

例11-9假设一信号采样频率为2000Hz，设计一个10阶低通ChebyshevI型数字滤波器，通带波纹为0.5dB，截止频率为300Hz。fs=2000;n=10;Rp=0.5;W1=300;Wn=W1/(fs/2);[b,a]=cheby1(n,Rp,Wn);%频率响应sampleN=512;freqz(b,a,sampleN,fs);axis([0500-300100]);11.3IIR滤波器《MATLAB程序设计》2611.3.3切比雪夫滤波器

cheby2函数cheby2函数用于设计切比雪夫数字滤波器中chebyshevII型。其调用格式为：[b，a]=cheby2(n，Rs，Ws)[b，a]=cheby2(n，Rs，Ws，ftype)[z，p，k]=cheby2(___)[A，B，C，D]=cheby2(___)[___]=cheby2(___，'s')11.3IIR滤波器《MATLAB程序设计》2711.3.3切比雪夫滤波器

例11-10假设一信号采样频率为2000Hz，设计一个10阶低通ChebyshevⅡ型数字滤波器，阻带衰减为20dB，通带波纹为0.5dB，截止频率为300Hz。fs=2000;n=10;Rs=20;W1=300;Wn=W1/(fs/2);[b,a]=cheby2(n,Rs,Wn);%频率响应sampleN=512;freqz(b,a,sampleN,fs);axis([0500-8020]);11.3IIR滤波器《MATLAB程序设计》2811.3.4椭圆滤波器

ellip函数用于设计椭圆数字滤波器。其调用格式为：[b，a]=ellip(n，Rp，Rs，Wp)[b，a]=ellip(n，Rp，Rs，Wp，ftype)[z，p，k]=ellip(___)[A，B，C，D]=ellip(___)[___]=ellip(___，'s')11.3IIR滤波器《MATLAB程序设计》2911.3.4椭圆滤波器

例11-11假设一信号采样频率为2kHz，设计一个8阶低通椭圆数字滤波器，阻带衰减50dB，通带波纹3dB，截止频率为300Hz。fs=2000;n=8;Rp=3;Rs=50;W1=300;Wn=W1/(fs/2);[b,a]=ellip(n,Rp,Rs,Wn);%频率响应sampleN=512;freqz(b,a,sampleN,fs);axis([0500-8020]);11.3IIR滤波器《MATLAB程序设计》3011.3.5信号分析实例

例11-12假设一个正弦信号包含50Hz和150Hz两个振动频率成分，采用IIR滤波器对该正弦信号进行处理，分别设计高通和低通滤波器，分别提取其低于100Hz和高于100Hz的信号成分。11.4FIR滤波器《MATLAB程序设计》3111.4.1FIR滤波器基本原理

FIR滤波器系统函数为，其中h(n)表示系统脉冲响应函数。11.4FIR滤波器《MATLAB程序设计》3211.4.2基于窗函数设计FIR滤波器

基于窗函数的FIR数字滤波器设计方法较为简单，其主要步骤如下：（1）由数字滤波器理想特性进行傅里叶逆变换，获得理想滤波器单位脉冲响应。一般假设理想低通滤波器的截止频率为ωc，其幅频特性为（2）由性能指标确定窗函数W(n)和窗口长度N。（3）求滤波器的单位脉冲响应h(n)（4）检验滤波器性能指标。11.4FIR滤波器《MATLAB程序设计》3311.4.2基于窗函数设计FIR滤波器

例11-13用窗函数设计一个线性相位FIR低通滤波器，满足性能指标为：通带边界频率wp=0.5π，阻带频率ωs=0.76π。Nw=N;wc=(wp+ws)/2;n=0:N-1;alpha=(N-1)/2;m=n-alpha+0.00001;hd=sin(wc*m)./(pi*m);win=hanning(Nw);h=hd.*win';b=h;freqz(b,1,512)wp=0.5*pi;ws=0.76*pi;width=ws-wp;N=ceil(8*pi/width);if(rem(N,2))==0N=N+1;end11.4FIR滤波器《MATLAB程序设计》3411.4.2基于窗函数设计FIR滤波器

例11-14用fir1设计一个32阶的FIR带通滤波器，通带为0.35≤ω≤0.65。n=32;W1=0.35;W2=0.65;Wn=[W1,W2];b=fir1(n,Wn);%频率响应sampleN=512;freqz(b,1,sampleN);axis([01-10050])11.4FIR滤波器《MATLAB程序设计》3511.4.2基于窗函数设计FIR滤波器

例11-15设计一个32阶的高通滤波器，截止频率为0.48，窗函数选用30dB波纹的chebwin窗函数。n=32;W1=0.48;nw=n+1;Rs=30;b=fir1(n，W1，'high'，chebwin(nw，Rs));%频率响应sampleN=512;freqz(b，1，sampleN);11.4FIR滤波器《MATLAB程序设计》3611.4.2基于窗函数设计FIR滤波器

例11-16设计一个32阶低通滤波器并绘制其期望频率响应与实际频率响应曲线。n=32;f=[00.60.61];m=[1100];b=fir2(n，f，m);%频率响应sampleN=512;[h，w]=freqz(b，1，sampleN);plot(f，m，'-');holdonplot(w/pi，abs(h)，'-.')xlabel('归一化频率');ylabel('幅值')ylim([-0.1，1.2])legend('期望频率响应'，'实际频率响应')11.4FIR滤波器《MATLAB程序设计》3711.4.3FIR滤波器优化设计

MATLAB信号处理工具箱提供了比基于窗函数法FIR滤波器设计函数fir1、fir2更为通用的函数firls和remez。这两个函数包采用不同的优化方法进行FIR滤波器设计。firls是fir1和fir2的扩展，基本原则是利用最小二乘法使期望的频率响应和实际频率响应间的误差最小。函数remez则通过采用Parks-McCellan交换算法和Chebyshev近似理论来设计FIR滤波器，使实际频率响应达到最优拟合期望频率响应的目的。11.4FIR滤波器《MATLAB程序设计》3811.4.3FIR滤波器优化设计

f=[00.30.40.60.70.9];a=[01000.50.5];n=32;%采用firls函数包进行设计b=firls(n，f，a);sampleN=512;[h，w]=freqz(b，1，sampleN);figureplot(f，a，'-');holdon例11-17分别用firls和remez函数包设计一个32阶的段线性带通滤波器，其理想幅频响应为f=[0，0.3，0.4，0.6，0.7，0.9]，a=[0，1，0，0，0.5，0.5]，要求完成滤波器设计并绘制期望的实际频率响应。plot(w/pi，abs(h)，'-.');xlabel('归一化频率');ylabel('幅值')legend('期望频率响应'，'实际频率响应')%采用remez函数包进行设计b2=remez(n，f，a);[h2，w2]=freqz(b，1，sampleN);figureplot(f，a，'-');holdonplot(w2/pi，abs(h2)，'-.');xlabel('归一化频率');ylabel('幅值')legend('期望频率响应'，'实际频率响应')11.4FIR滤波器《MATLAB程序设计》3911.4.4信号分析实例

例11-18本节将以MATLAB自带的chirp信号分析为例，分别设计FIR低通、高通滤波器，完成chirp中高频和低频成分提取。已知chirp信号采样频率fs为8192Hz，其主要频率成分集中在fs/4频率以上。分别设计一个阶数为34的低通和高通滤波器，以fs/4为频率界限，分别实现chirp低于fs/4和高于fs/4频率处成分提取。11.5特殊波形发生函数《MATLAB程序设计》4011.5.1扫频信号

MATLAB中的chirp函数用于产生扫频余弦信号，其调用格式如下：y=chirp(t，f0，t1，f1)y=chirp(t，f0，t1，f1，method)y=chirp(t，f0，t1，f1，method，phi)y=chirp(t，f0，t1，f1，'quadratic'，phi，shape)y=chirp(___，cplx)11.5特殊波形发生函数《MATLAB程序设计》4111.5.1扫频信号

例11-19模拟仿真一个线性扫频信号，其采样时长为1s，采样频率为1kHz，要求其初始时刻频率值为0Hz，在0.5s时，其频率值为250Hz。%参数设置fs=1e3;T=1/fs;f1=0;f2=50;t1=0;t2=1;t=t1:T:t2;tx=0.5;%模拟仿真y=chirp(t，f1，tx，f2);%信号显示figureplot(t，y);xlabel('t');ylabel('y');11.5特殊波形发生函数《MATLAB程序设计》4211.5.1扫频信号

例11-19模拟仿真一个线性扫频信号，其采样时长为1s，采样频率为1kHz，要求其初始时刻频率值为0Hz，在0.5s时，其频率值为250Hz。figurepspectrum(y，fs，'spectrogram')11.5特殊波形发生函数《MATLAB程序设计》4311.5.1扫频信号

例11-20模拟仿真一个二次扫频信号，其采样时长为1s，采样频率为1kHz，要求其初始时刻频率值为0Hz，在0.5s时，其频率值为250Hz。%参数设置fs=1e3;T=1/fs;f1=0;f2=250;t1=0;t2=1;t=t1:T:t2;tx=0.5;%模拟仿真y=chirp(t，f1，tx，f2，'quadratic');%信号显示figureplot(t，y);xlabel('t');ylabel('y');figurepspectrum(y，fs，'spectrogram')11.5特殊波形发生函数《MATLAB程序设计》4411.5.2冲击信号

冲击函数的数学定义如下，其中N为非零整数。11.5特殊波形发生函数《MATLAB程序设计》4511.5.2冲击信号

例11-21模拟仿真冲击函数，令N取值为13，