信号处理仿真实验
目录
一、离散傅立叶变换 (3)
1、补零DFT (3)
2、高密度频谱和高分辨率频谱 (4)
3、用DFT进行谱分析 (6)
二、数字滤波器的设计 (11)
4、IIR滤波器的设计 (11)
5、FIR滤波器的设计 (14)
三、短时傅立叶变换分析信号特性 (17)
6、短时傅立叶变换分析信号特性 (17)
四、功率谱估计与应用 (18)
7、采用周期图法估计信号的功率谱 (18)
8、采用平均周期图法估计信号的功率谱 (20)
9、采用分段平均周期图法估计信号的功率谱 (21)
五、同态信号处理 (23)
10、同态信号处理 (24)
六、非平稳信号分析与处理 (26)
七、Kalman滤波在机动目标跟踪中的应用 (31)
八、机载R-D成像................................................................................ 错误！未定义书签。
一、 离散傅立叶变换
1、补零DFT
5()()x n R n =，求N 0分别取10，20时的X(k)。
解：设N 0=10，则
00
22214
10
5
sin()
2()()0,1,...,9sin()10N j kn j kn j k
N n n k X k x n e
e
e
k k πππππ
----===
===∑∑
设N 0=20，则
00
2214
20
5
sin()
4
()()0,1,...,19sin(
)
20
N j
kn j kn j k N n n k X k x n e
e
e
k k πππ
ππ
----===
===∑∑
用Matlab 进行仿真，得图1-1。
Matlab 源程序：
%X(k)与X(ejw)的关系
clear;
n=0:4;x=[ones(1,5)]; k=0:999;w=(pi/500)*k;
X=x*(exp(-j*pi/500)).^(n'*k);%计算离散时间傅立叶变换X(ejw) Xe=abs(X);
subplot(321);stem(n,x);ylabel('x(n)');
subplot(322);plot(w/pi,Xe);ylabel('|X(ejw)|'); N=15;x=[ones(1,5),zeros(1,N-5)]; n=0:N-1; k=0:9;
X=x*(exp(-j*pi/5)).^(n'*k);%N=10点离散傅立叶变换 magX=abs(X);
subplot(323);stem(n,x);ylabel('x(n)');
subplot(324);stem(k,magX);axis([0,10,0,5]);ylabel('|X(k)|'); N=20;x=[ones(1,5),zeros(1,N-5)]; n=0:1:N-1; k=0:N-1;
X=x*(exp(-j*pi/10)).^(n'*k);;%N=20点离散傅立叶变换 magX=abs(X);
subplot(325);stem(n,x);ylabel('x(n)');
subplot(326);stem(k,magX);axis([0,20,0,5]);ylabel('|X(k)|');
图1-1 ()X k 与()j X e ω
的关系
结论：
（1） 填零是给原始序列的运算，给给原始序列的离散时间傅立叶变换提供间隔
较密的样本。
（2） 为画出()j X e ω
，只需要5点的()X k 用内插公式即可得到()j X e ω
，但实
际上是用10或20点的()X k 来填充()j X e ω
的值。
（3） 填零运算提供了较密的频谱，而没有增加任何新的信息，因此它不能提供
高分辨率的频谱。
（4） 为得到高分辨率的频谱，需从实验或观察中取得更多的数据。
2、高密度频谱和高分辨率频谱
为了说明高密度频谱和高分辨率频谱之间的区别，考察序列
()2cos(0.35)cos(0.5)x n n n ππ=+
（1） 当010n ≤≤时，确定并画出x(n)的离散傅立叶变化。
（2） 当()
010()0
1040
x n n x n ≤
≤?时，确定并画出x(n)的离散傅立叶变换。
（3） 当040n ≤≤时，确定并画出x(n)的离散傅立叶变换。
用Matlab 进行仿真，得图1-2。 Matlab 源程序：
%高密度频谱与高分辨率频谱
clear;
figure(1);
N=40;M=10;
n=0:N-1;
x=2*cos(0.35*pi*n)+cos(0.5*pi*n);
n1=0:M-1;x1=x(1:M);
subplot(231);stem(n1,x1);title('没有足够采样点的信号');
k1=0:M-1;w1=2*pi*k1/M;
X1=x1*(exp(-j*2*pi/M)).^(n1'*k1);%计算M点DFT
Y1=abs(X1);
subplot(234);stem(w1/pi,Y1);title('信号的频谱');
n2=0:N-1;x2=[x1,zeros(1,N-M)];
subplot(232);stem(n2,x2);title('填零信号');
k2=0:N-1;w2=2*pi*k2/N;
X2=x2*(exp(-j*2*pi/N)).^(n2'*k2);%计算N点DFT
Y2=abs(X2);
subplot(235);stem(w2/pi,Y2);title('高密度频谱');
subplot(233);stem(n,x);title('有足够采样点的信号');
k=0:N-1;w=2*pi*k/N;
X=x*(exp(-j*2*pi/N)).^(n'*k);%计算N点DFT
Y=abs(X);
subplot(236);stem(w/pi,Y);title('高分辨率频谱');
图1-2 高密度频谱与高分辨率频谱结论：
（1） 当010n ≤≤时的序列x(n)与X(k)，从X(k)图中几乎无法看出有关信号的
频谱的信息。
（2） 将x(n)补30个零点时的x(n)与X(k)，这时的谱线相当密，但从中很难看出
信号的频率成分，故称作高密度频谱。
将x(n)的长度加长到40时的x(n)与X(k)，这时可以清晰的看出信号的频谱成分（120.35,0.5ωπωπ==），故称为高分辨率频谱。
3、用DFT 进行谱分析
信号的谱分析就是计算信号的傅立叶变换。连续信号与系统的傅立叶变换不便于直接用计算机进行计算，其应用受到限制。DFT 是离散信号与系统分析的有力工具，对连续信号与系统可以通过时域采样，应用DFT 进行近似谱分析。
对模拟信号()2sin(4)5cos(8)
a x t t t ππ=+以间隔T 进行采样，得到N 点序列()x n ，用N 点DFT 得到()a x t 幅度谱的估计。
（1） T=0.01s ，N=40或N=50，一个能提供精确的()a x t 的幅度谱，画出DFT 的幅度谱。 （2） T=0.005s ，N=40或N=50，画出DFT 的幅度谱。
用Matlab 进行仿真，得图1-3。 Matlab 源程序：
%用DFT 进行频谱分析 T=0.01;N=40;n=0:N-1;t=n*T; xn=2*sin(4*pi*t)+5*cos(8*pi*t); Xk=fft(xn,N); magXk=abs(Xk); k=0:N-1;
subplot(241);plot(t,xn);axis([0,0.4,-7.5,7]); title('T=0.01s,t=0.4s');ylabel('x(t)');
subplot(245);stem(k,magXk);title('T=0.01s,N=40');ylabel('X(k)');
T=0.01;N=50;n=0:N-1;t=n*T; xn=2*sin(4*pi*t)+5*cos(8*pi*t); Xk=fft(xn,N); magXk=abs(Xk); k=0:N-1;
subplot(242);plot(t,xn);axis([0,0.5,-7.5,7]); title('T=0.01s,t=0.5s');
subplot(246);stem(k,magXk);title('T=0.01s,N=50');
T=0.005;N=40;n=0:N-1;t=n*T; xn=2*sin(4*pi*t)+5*cos(8*pi*t); Xk=fft(xn,N);
magXk=abs(Xk); k=0:N-1;
subplot(243);plot(t,xn);axis([0,0.2,-7.5,7]); title('T=0.005s,t=0.2s');
subplot(247);stem(k,magXk);title('T=0.005s,N=40');
T=0.005;N=50;n=0:N-1;t=n*T; xn=2*sin(4*pi*t)+5*cos(8*pi*t); Xk=fft(xn,N); magXk=abs(Xk); k=0:N-1;
subplot(244);plot(t,xn);axis([0,0.25,-7.5,7]); title('T=0.005s,t=0.25s');
subplot(248);stem(k,magXk);title('T=0.005s,N=50');
图1-3 用DFT 进行频谱分析
结论：采样间隔T=0.005s ，采样点数N=50时的幅度谱是最精确的()a x t 的幅度谱，而其他情况都存在频谱泄漏。
用DFT 分析连续时间信号的频谱。已知一模拟信号()()t
a x t e u t -=，现以采样率
20s f Hz =进行采样。用DFT 计算当序列长度L=100和L=20时，N=200点的幅度频
谱样值并通过作图与理论上准确的频谱样值进行比较。 解：原信号的傅立叶变换
1
()()1j t a a X j x t e dt j ∞
-Ω-∞
Ω=
=
+Ω
?
其幅度为
|()|a X j Ω=
1-4所示。用DFT 计算的频谱如图1-5所示。
Matlab 源程序为：
%理论上的频谱 figure(1);
Omeger=0:0.1:20*pi; Xa=1./sqrt(1+Omeger.^2);
subplot(111);plot(Omeger/pi,Xa);title('|Xa(j*Omeger)|'); %用DFT 计算的频谱 fs=20;
L=100;N=200;n=0:L-1;t1=n/fs;
xn1=exp(-t1);xn=[xn1,zeros(1,N-L)]; Xk1=fft(xn,N);magXk1=abs(Xk1); k1=0:N-1;
L=20;N=200;n=0:L-1;t2=n/fs;
xn2=exp(-t2);xn=[xn2,zeros(1,N-L)]; Xk2=fft(xn,N);magXk2=abs(Xk2); k2=0:N-1; figure(2);
subplot(221);plot(t1,xn1);title('xa(t) t=5s'); subplot(223);plot(t2,xn2);title('xa(t) t=1s');
subplot(222);plot(k1,magXk1);title('X(k) L=100 N=200'); subplot(224);plot(k2,magXk2);title('X(k) L=20 N=200');
图1-4 理论上的频谱 图1-5 用DFT 计算的频谱
结论：
（1） 当序列长度为100时，进行200点DFT 计算的结果混叠与泄漏影比较小，
基本上接近原信号的频谱。因为按给定的20s f Hz =，相当于取信号的最高频率20h f Hz =，故在[0,]h f 频率范围内的信号能量为
202
201()0.4952h E X j d π
π
π
-=
ΩΩ=?
而信号总的能量为
2
11
()arctan 0.522x E X j d π
π∞
∞
-∞
-∞
=
ΩΩ=Ω=?
所以/99%h x E E =，基本上满足不混叠的要求。
（2） 当序列长度为20时，进行200点DFT 计算，由于截取x(n)的长度太短，
1
()0.30790LT t LT x t e e
-===
=>> 所以频谱因泄漏出现较大的波动，以致与原信号频谱有较大的差别。
用加窗函数DFT 分析离散时间信号的频谱。已知信号为0.05()n
x n e
-=，分别加矩形窗、
汉宁窗和三角窗对其进行DFT 。用Matlab 画得如图1-6所示结果。 Matlab 源程序如下：
%加窗后的信号及频谱 figure(1);
L=20;N=200;n=0:L-1;
xn=exp(-0.05*n);xn1=[xn,zeros(1,N-L)];
Xk=fft(xn1,N);magXk=abs(Xk);
k=0:N-1;
subplot(211);stem(k,xn1);title('矩形窗x(n)'); subplot(212);stem(k,magXk);title('矩形窗X(k)'); figure(2);
w=(hanning(L))';xn1=xn.*w;
xn2=[xn1,zeros(1,N-L)];
subplot(211);stem(k,xn2);title('汉宁窗x(n)');
Xk=fft(xn2,N);magXk=abs(Xk);
subplot(212);stem(k,magXk);title('汉宁窗X(k)'); figure(3);
w=(triang(L))';xn1=xn.*w;
xn2=[xn1,zeros(1,N-L)];
subplot(211);stem(k,xn2);title('三角窗x(n)');
Xk=fft(xn2,N);magXk=abs(Xk);
subplot(212);stem(k,magXk);title('汉宁窗X(k)'); figure(4);
w=(hamming(L))';xn1=xn.*w;
xn2=[xn1,zeros(1,N-L)];
subplot(211);stem(k,xn2);title('汉明窗x(n)');
Xk=fft(xn2,N);magXk=abs(Xk);
subplot(212);stem(k,magXk);title('汉明窗X(k)'); figure(5);
w=(blackman(L))';xn1=xn.*w;
xn2=[xn1,zeros(1,N-L)];
subplot(211);stem(k,xn2);title('布拉克曼x(n)');
Xk=fft(xn2,N);magXk=abs(Xk);
subplot(212);stem(k,magXk);title('布拉克曼X(k)');
图1-6 加窗后的信号及频谱
结论：
（1）实际工作中往往要把信号的观察时间限制在一定的时间间隔内，就需要用有限个数据，即将信号截断，相当于将信号乘以窗函数。在时域的截断，相当于所
研究的频谱与窗函数频谱的周期卷积，周期卷积的结果造成频谱从其正常频谱
扩展开来，成为“泄漏”。泄漏效应是DFT固有的，可用加窗函数进行抑制。用
加窗函数加权使有限长度的输入信号周期延拓后，在边界上尽量减少不连续程
度，从而达到抑制频谱泄漏的目的。
（2）从图1-6可以看出，加不同窗函数得到的过渡带由窄到宽的顺序为：
矩形窗→三角窗→汉宁窗→汉明窗→布拉克曼窗
加不同窗函数得到的最小阻带衰减由低到高的顺序为：
矩形窗→三角窗→汉宁窗→汉明窗→布拉克曼窗
可见，若要得到最窄的过渡带，选择矩形窗，要得到最高的最小阻带衰减，选布拉克曼窗，而若要兼顾过渡带和最小阻带衰减，则选择三角窗、汉宁窗和
汉明窗。
二、数字滤波器的设计
4、IIR滤波器的设计
用并联型实现系统函数为：
()4
3214
32142121618271131----------+++++-=
z z z z z z z z z H 的IIR 数字滤波器，求单位脉冲响应和单位阶跃信号的输出。 程序如下：
%并联型单位脉冲响应和单位阶跃信号的输出 b=[1,-3,11,27,18]; a=[16,12,2,-4,-1];
[C,B,A]=dir2par(b,a);%直接型转换成并联型 N=25;
delta=impseq(0,0,N);
h=parfiltr(C,B,A,delta);%并联型单位脉冲响应 x=[ones(1,5),zeros(1,N-5)];
y=parfiltr(C,B,A,x);%并联型单位阶跃信号的输出 subplot(121); stem(h);
title('并联型h(n)'); subplot(122); stem(y);
title('并联型y(n)');
function [C,B,A]=dir2par(b,a) %直接型转换成并联型 M=length(b); N=length(a);
[r1,p1,C]=residuez(b,a);%展开成部分分式
p=cplxpair(p1,10000000*eps);%寻找共轭复数对
I=cplxcomp(p1,p); %记下p1与p 的对应关系
r=r1(I); %使r 与p 对应（就像r1与p1对应一样） K=floor(N/2); B=zeros(K,2); A=zeros(K,3); if (K*2==N)
for i=1:2:N-2
Brow=r(i:i+1,:); Arow=p(i:i+1,:);
[Brow,Arow]=residuez(Brow,Arow,[]);%把共轭复数对对应的分式合并到一起 B(fix((i+1)/2),:)=real(Brow); A(fix((i+1)/2),:)=real(Arow); end
[Brow,Arow]=residuez(r(N-1),p(N-1),[]); B(K,:)=[real(Brow),0];A(K,:)=[real(Arow),0];
else
for i=1:2:N-1
Brow=r(i:i+1,:);
Arow=p(i:i+1,:);
[Brow,Arow]=residuez(Brow,Arow,[]);%把共轭复数对对应的分式合并到一起
B(fix((i+1)/2),:)=real(Brow);
A(fix((i+1)/2),:)=real(Arow);
end
end
function I=cplxcomp(p1,p)
N=length(p);
for k=1:N
I(k)=find(abs(p1-p(k))<10000000*eps);
End
function [x,n]=impseq(n0,n1,n2)
%产生单位脉冲
n=[n1:n2];
x=[(n-n0)==0];
function y=parfiltr(C,B,A,x)
%滤波器的并联型实现
[K,L]=size(B);
N=length(x);
w=zeros(K+1,N);
if (isempty(C)~=1) %C可能是空矩阵
w(1,:)=filter(C,1,x);
end
for i=1:K
w(i+1,:)=filter(B(i,:),A(i,:),x);%各个分支分别滤波
end
y=sum(w); %各个分支滤波结果求和
单位脉冲响应和单位阶跃信号的输出
结果：
5、FIR 滤波器的设计
已知数字滤波器的性能指标如下所示：
0.20.250.350p s s Rp dB A dB
ωπωπ
====
用窗函数法设计一个FIR 低通滤波器。 选择汉明窗。 程序如下：
%用Hamming 窗设计低通滤波器 wp=0.2*pi; ws=0.3*pi;
tr_width=ws-wp;
M=ceil(6.6*pi/tr_width)+1;
disp(['滤波器的长度为',num2str(M)]); n=0:M-1;
%理想LPF 的截止频率 wc=(ws+wp)/2; hd=ideal_lp(wc,M); w_ham=(hamming(M))'; h=hd.*w_ham;
[db,mag,pha,grd,w]=freqz_m(h,[1]); delta_w=2*pi/1000;
%求出实际通带波动
Rp=-(min(db(1:1:wp/delta_w+1)));
disp(['实际通带波动为',num2str(Rp)]);
%求出最小阻带衰减
As=-round(max(db(ws/delta_w+1:1:501))); disp(['最小阻带衰减为',num2str(As)]);
subplot(111)
subplot(221)
stem(n,hd);
title('理想冲激响应');
axis([0 M-1 -0.1 0.3]);
ylabel('hd(n)');
subplot(222)
stem(n,w_ham);
title('汉明窗');
axis([0 M-1 0 1.1]);
ylabel('w(n)');
subplot(223)
stem(n,h);
title('实际冲激响应');
axis([0 M-1 -0.1 0.3]);
xlabel('n');
ylabel('h(n)');
subplot(224)
plot(w/pi,db);
title('幅度响应(dB)');
axis([0 1 -100 10]);
grid;
xlabel('以pi为单位的频率');
ylabel('分贝数');
function hd=ideal_lp(wc,M)
%理想低通滤波器计算
%
%hd为0到M-1之间的理想脉冲响应
%wc为截止频率
%M为理想滤波器的长度
%
alpha=(M-1)/2;
n=0:M-1;
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
function [db,mag,pha,grd,w]=freqz_m(b,a)
%求取系统响应
%
%db为相对振幅
%mag为绝对振幅
%pha为相位响应
%grd为群延时
%w为频率样本点向量
%b为Ha(z)分子多项式系数（对FIR而言，b=h）%a为Hz(z)分母多项式系数（对FIR而言，a=1）%
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:501))';
w=(w(1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);
结果：
滤波器的长度为67
实际通带波动为0.03936
最小阻带衰减为52
三、 短时傅立叶变换分析信号特性
6、短时傅立叶变换分析信号特性
数字信号300,1,...,0),(=-=N N n n x ，在（50，100）和（200，260）分别有两个频率不同的正弦信号，用短时傅立叶变换分析其特性。
附程序:
function stft N=300;t=0:N-1; x=zeros(1,N);
x(50:100)=cos(pi*(t(50:100)-50)/10); x(200:260)=cos(pi*(t(200:260)-200)/5); subplot(2,2,1); plot(x); title('x'); X=fft(x); X=fftshift(X); subplot(2,2,2); plot(abs(X)); title('abs(X)'); Nw=20;
L=Nw/2;Ts=(N-Nw)/L+1; nfft=32;
TF=zeros(Ts,nfft); for i=1:Ts
xw=x((i-1)*L+1:i*L+L); temp=fft(xw,nfft); temp=fftshift(temp); TF(i,:)=temp; end
subplot(2,2,3);mesh(abs(TF)); title('abs(TF)');
subplot(2,2,4);contour(abs(TF)); title('abs(TF)');
实验结果
图 3－1
实验分析
由实验结果可知，在信号的频谱中可以看出观测到的信号中存在两种频率，但是不能得出它们的时间持续范围，而从信号的时频平面中不仅可以看出这两种信号的频率，而且还能够得到它们的时间持续范围。不确定性原理表明：对于短时傅立叶变换，其时间分辨力和频率分辨力是互相矛盾的，信号的时域波形与频谱不能同时获得高的分辨力。在实际信号分析中，应采取恰当的时间窗函数，由具体问题本身的特性来决定。
四、 功率谱估计与应用
这里这讨论周期图法。
周期图法是直接将离散信号x(n)进行傅立叶变换来求取功率谱估计的，即
2
1()[()],0,1,...,1S k DFT x n k N N
=
=-
由于DFT[x(n)]的周期为N ，因而求得的功率谱估计也是以N 为周期的，故称这种方法为周期图法。周期图法之所以广泛使用，是由于它与序列的频谱有对应关系，可以采用FFT 算法来快速计算。
7、采用周期图法估计信号的功率谱
求信号
12()sin(2)2cos(2)(),1,...,x n f n f n w n n N ππ=++=的功率谱。其中，
120.05,0.120,f f ==()n ω为白噪声。
（1） 信号长度N=256；
（2）信号长度N=1024。
Matlab程序如下：
%采用周期图法估计信号的功率谱
%
clf;
%Case 1:N=256;
N=256;
n=[0:N-1];
f1=0.05;
f2=0.120;
wn=randn(1,N);
xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn;
Pxx=10*log10(abs(fft(xn).^2)/(N+1));
f=(0:length(Pxx)-1)/length(Pxx);
figure(1);
plot(f,Pxx);
xlabel('Frequency');
ylabel('Power Spectrum(dB)');
title('Periodogram N=256');
grid;
%======================================================
%Case 2:N=1024;
N1=1024;
n=[0:N1-1];
f1=0.05;
f2=0.120;
wn1=randn(1,N1);
xn1=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn1;
Pxx1=10*log10(abs(fft(xn1).^2)/(N1+1));
f=(0:length(Pxx1)-1)/length(Pxx1);
figure(2);
plot(f,Pxx1);
xlabel('Frequency');
ylabel('Power Spectrum(dB)');
title('Periodogram N=1024');grid;
运行得如图4－1所示的功率谱估计。
结论：由图4－1可以看出，在频率0.05和0.12处功率谱有两个峰值，说明信号中含有0.05和0.12的周期成分。但功率谱都在很大范围内波动，而且并没有因信号取样点数N的增加而有明显改进。我们已经知道，N的增加并不会使周期图的方差趋于零，这样N越大，周期图上的随机起伏就更加密集和剧烈。
功率谱估计中，有限长的数据可视为无限长随机序列经矩形窗截断的结果，这也意味着对自相关函数的加窗，使得功率谱与窗函数卷积，这种频域卷积会产生频谱泄漏，容易使弱信号的主瓣被强信号所淹没，造成频谱的模糊和失真，使得周期图分辨率较低。
（a ） (b)
图4－1 周期图谱估计
用有限长样本序列的FFT 来表示随机序列的功率谱是一种估计或近似，不可避免存在着误差。为了减少误差，使功率谱估计更加平滑，可采样两种改进的周期图法：分段平均周期图法和加窗平均周期图法。
8、采用平均周期图法估计信号的功率谱
求信号
12()sin(2)2cos(2)(),1,...,x n f n f n w n n N ππ=++=的功率谱。其中，
120.05,0.120,f f ==()n ω为白噪声，信号程度N=1024。用两种平均周期图法求信号
的功率谱估计。 Matlab 程序如下：
%采用平均周期图法估计信号的功率谱 % clf;
N=1024;Ns=256;n=[0:N-1];
f1=0.05;f2=0.120;wn=randn(1,N);
xn=sin(2*pi*f1*n)+2*cos(2*pi*f2*n)+wn; pxx1=abs(fft(xn(1:256),Ns).^2)/Ns; pxx2=abs(fft(xn(257:512),Ns).^2)/Ns; pxx3=abs(fft(xn(513:768),Ns).^2)/Ns; pxx4=abs(fft(xn(769:1024),Ns).^2)/Ns;
Pxx=10*log10((pxx1+pxx2+pxx3+pxx4)/4); f=(0:length(Pxx)-1)/length(Pxx); figure(1);plot(f,Pxx); xlabel('Frequency');
ylabel('Power Spectrum(dB)');
title('Averaged Periodogram(no overlap) N=4*256');grid; %Estimate PSD by the halfoverlapping signalsection
%====================================================== N=1024;Ns=256;n=[0:N-1];
f1=0.05;f2=0.120;wn=randn(1,N);
目录 摘要 (2) 第一章绪论 (3) 1.1 语音课设的意义 (3) 1.2 语音课设的目的与要求 (3) 1.3 语音课设的基本步骤 (3) 第二章设计方案论证 (5) 2.1 设计理论依据 (5) 2.1.1 采样定理 (5) 2.1.2 采样频率 (5) 2.1.3 采样位数与采样频率 (5) 2.2 语音信号的分析及处理方法 (6) 2.2.1 语音的录入与打开 (6) 2.2.2 时域信号的FFT分析 (6) 2.2.3 数字滤波器设计原理 (7) 2.2.4 数字滤波器的设计步骤 (7) 2.2.5 IIR滤波器与FIR滤波器的性能比较 (7) 第三章图形用户界面设计 (8) 3.1 图形用户界面概念 (8) 3.2 图形用户界面设计 (8) 3.3 图形用户界面模块调试 (9) 3.3.1 语音信号的读入与打开 (9) 3.3.2 语音信号的定点分析 (9) 3.3.3 N阶高通滤波器 (11) 3.3.4 N阶低通滤波器 (12) 3.3.5 2N阶带通滤波器 (13) 3.3.6 2N阶带阻滤波器 (14) 3.4 图形用户界面制作 (15) 第四章总结 (18) 附录 (19) 参考文献 (24)
摘要 数字信号处理是将信号以数字方式表示并处理的理论和技术。数字信号处理与模拟信号处理是信号处理的子集。 数字信号处理的目的是对真实世界的连续模拟信号进行测量或滤波。因此在进行数字信号处理之前需要将信号从模拟域转换到数字域，这通常通过模数转换器实现。而数字信号处理的输出经常也要变换到模拟域，这是通过数模转换器实现的。 数字信号处理的算法需要利用计算机或专用处理设备如数字信号处理器（DSP）和专用集成电路（ASIC）等。数字信号处理技术及设备具有灵活、精确、抗干扰强、设备尺寸小、造价低、速度快等突出优点，这些都是模拟信号处理技术与设备所无法比拟的。 数字信号处理的核心算法是离散傅立叶变换(DFT)，是DFT使信号在数字域和频域都实现了离散化，从而可以用通用计算机处理离散信号。而使数字信号处理从理论走向实用的是快速傅立叶变换(FFT)，FFT的出现大大减少了DFT的运算量，使实时的数字信号处理成为可能、极大促进了该学科的发展。 MATLAB是矩阵实验室（Matrix Laboratory）的简称，和Mathematica、Maple 并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完相同的事情简捷得多，并且mathwork也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++ ，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。
做脑电波信号处理滴嘿嘿。。Matlab addicted Codes %FEATURE EXTRACTER function [features] = EEGfeaturetrainmod(filename,m) a = 4; b = 7; d = 12; e = 30; signals = 0; for index = 1:9; % read in the first ten EEG data because the files are numbered as ha11test01 rather than ha11test1. s = [filename '0' num2str(index) '.dat']; signal = tread_wfdb(s); if signals == 0; signals = signal; else signals = [signals signal]; end end for index = 10:1:m/2; % read in the rest of the EEG training data s = [filename num2str(index) '.dat']; signal = tread_wfdb(s); if signals == 0;
signals = signal; else signals = [signals signal]; end end %%%%% modification just for varying the training testing ratio ------ for index = 25:1:25+m/2; % read in the rest of the EEG training data s = [filename num2str(index) '.dat']; signal = tread_wfdb(s); if signals == 0; signals = signal; else signals = [signals signal]; end end %%%%%end of modification just for varying the training testing ratio----- for l = 1:m % exrating features (power of each kind of EEG wave forms) [Pxx,f]=pwelch(signals(:,l)-mean(signals(:,l)), [], [], [], 200); % relative power fdelta(l) = sum(Pxx(find(f
安装好MATLAB 2012后再安装目录下点击setup.exe 会出现 "查找安装程序类时出错，查找类时出现异常"的错误提示。该错误的解决方法是进入安装目录下的bin 文件夹双击matlab.exe 对安装程序进行激活。这是可以对该matlab.exe 创建桌面快捷方式，以后运行程序是直接双击该快捷方式即可。 信号运算 1、 信号加 MATLAB 实现: x=x1+x2 2、 信号延迟 y(n)=x(n-k) 3、 信号乘 x=x1.*x2 4、 信号变化幅度 y=k*x 5、 信号翻转 y=fliplr(x) 6、 信号采样和 数学描述：y=∑=2 1)(n n n n x MATLAB 实现： y=sum(x(n1:n2)) 7、 信号采样积 数学描述：∏==2 1)(n n n n x y MATLAB 实现： y=prod(x(n1:n2)) 8、 信号能量 数学描述：∑∞ -∞ == n x n x E 2 | )(| MATLAB 实现：Ex=sum(abs(x)^2)
9、 信号功率 数学描述：∑-== 1 2 | )(|1 P N n x n x N MATLAB 实现：Px=sum((abs(x)^2)/N MATLAB 窗函数 矩形窗 w=boxcar(n) 巴特利特窗 w=bartlett(n) 三角窗 w=triang(n) 布莱克曼窗 w=blackman(n) w=blackman(n,sflag) 海明窗 w=haiming(n) W=haiming(n,sflag) sflag 用来控制窗函数首尾的两个元素值，其取值为symmetric 、periodic 汉宁窗 w=hanning(n) 凯塞窗 w=Kaiser(n,beta) ,beta 用于控制旁瓣的高度。n 一定时，beta 越大，其频谱的旁瓣越小，但主瓣宽度相应增加；当beta 一定时，n 发生变化，其旁瓣高度不变。 切比雪夫窗：主瓣宽度最小，具有等波纹型，切比雪夫窗在边沿的采样点有尖峰。 W=chebwin(n,r)
数字信号处理Matlab 实现实例 第１章离散时间信号与系统 例1-1 用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。 解 MATLAB程序如下： a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 图1.1给出了卷积结果的图形，求得的结果存放在数组c中为：{-2 -4 1 3 1 5 1 -3}。 例1-2 用MATLAB计算差分方程 当输入序列为时的输出结果。 解 MATLAB程序如下： N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)];
k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n');ylabel('幅度') 图 1.2 给出了该差分方程的前41个样点的输出，即该系统的单位脉冲响应。 例1-3 用MATLAB 计算例1-2差分方程 所对应的系统函数的DTFT 。 解 例1-2差分方程所对应的系统函数为： 123 123 0.80.440.360.02()10.70.450.6z z z H z z z z -------++= +-- 其DTFT 为 23230.80.440.360.02()10.70.450.6j j j j j j j e e e H e e e e ωωωω ωωω--------++= +-- 用MATLAB 计算的程序如下： k=256; num=[0.8 -0.44 0.36 0.02]; den=[1 0.7 -0.45 -0.6]; w=0:pi/k:pi; h=freqz(num,den,w); subplot(2,2,1); plot(w/pi,real(h));grid title('实部') xlabel('\omega/\pi');ylabel('幅度')
系（院）物理与电子工程学院专业电子信息工程题目语音信号的处理与分析 学生姓名 指导教师 班级 学号 完成日期：2013 年5 月 目录 1 绪论.............................................................................................................. 错误!未定义书签。 1.1课题背景及意义................................................................................. 错误!未定义书签。 1.2国内外研究现状................................................................................. 错误!未定义书签。 1.3本课题的研究内容和方法................................................................. 错误!未定义书签。 1.3.1 研究内容................................................................................ 错误!未定义书签。 1.3.2 开发环境................................................................................ 错误!未定义书签。 2 语音信号处理的总体方案............................................................................ 错误!未定义书签。 2.1 系统基本概述.................................................................................... 错误!未定义书签。 2.2 系统基本要求与目的........................................................................ 错误!未定义书签。 2.3 系统框架及实现................................................................................ 错误!未定义书签。 2.3.1 语音信号的采样.................................................................... 错误!未定义书签。 2.3.2 语音信号的频谱分析............................................................ 错误!未定义书签。 2.3.3 音乐信号的抽取.................................................................... 错误!未定义书签。 2.3.4 音乐信号的AM调制.............................................................. 错误!未定义书签。 2.3.5 AM调制音乐信号的同步解调............................................... 错误!未定义书签。 2.4系统设计流程图................................................................................. 错误!未定义书签。 3 语音信号处理基本知识................................................................................ 错误!未定义书签。 3.1语音的录入与打开............................................................................. 错误!未定义书签。 3.2采样位数和采样频率......................................................................... 错误!未定义书签。 3.3时域信号的FFT分析......................................................................... 错误!未定义书签。 3.4切比雪夫滤波器................................................................................. 错误!未定义书签。 3.5数字滤波器设计原理......................................................................... 错误!未定义书签。 4 语音信号实例处理设计................................................................................ 错误!未定义书签。 4.1语音信号的采集................................................................................. 错误!未定义书签。
语音信号处理系统设计 摘要：语音信号处理是研究用数字信号处理技术对语音信号进行处理的一门学科。语音信号处理的目的是得到某些参数以便高效传输或存储,或者是用于某种应用，如人工合成出语音、辨识出讲话者、识别出讲话内容、进行语音增强等。本文简要介绍了语音信号采集与分析以及语音信号的特征、采集与分析方法，并在采集语音信号后，在MATLAB 软件平台上进行频谱分析,并对所采集的语音信号加入干扰噪声，对加入噪声的信号进行频谱分析，设计合适的滤波器滤除噪声，恢复原信号。利用MATLAB来读入（采集）语音信号，将它赋值给某一向量，再将该向量看作一个普通的信号，对其进行FFT变换实现频谱分析，再依据实际情况对它进行滤波，然后我们还可以通过sound命令来对语音信号进行回放，以便在听觉上来感受声音的变化。 关键词：Matlab，语音信号，傅里叶变换，滤波器 1课程设计的目的和意义 本设计课题主要研究语音信号初步分析的软件实现方法、滤波器的设计及应用。通过完成本课题的设计，拟主要达到以下几个目的： 1.1．了解Matlab软件的特点和使用方法。 1.2．掌握利用Matlab分析信号和系统的时域、频域特性的方法； 1.3．掌握数字滤波器的设计方法及应用。 1.4．了解语音信号的特性及分析方法。 1.5．通过本课题的设计，培养学生运用所学知识分析和解决实际问题的能力。 2 设计任务及技术指标 设计一个简单的语音信号分析系统，实现对语音信号时域波形显示、进行频谱分析，
利用滤波器滤除噪声、对语音信号的参数进行提取分析等功能。采用Matlab设计语言信号分析相关程序，并且利用GUI设计图形用户界面。具体任务是： 2.1．采集语音信号。 2.2．对原始语音信号加入干扰噪声，对原始语音信号及带噪语音信号进行时频域分析。 2.3．针对语音信号频谱及噪声频率，设计合适的数字滤波器滤除噪声。 2.4．对噪声滤除前后的语音进行时频域分析。 2.5.对语音信号进行重采样，回放并与原始信号进行比较。 2.6．对语音信号部分时域参数进行提取。 2.7．设计图形用户界面（包含以上功能）。 3 设计方案论证 3.1语音信号的采集 使用电脑的声卡设备采集一段语音信号，并将其保存在电脑中。 3.2语音信号的处理 语音信号的处理主要包括信号的提取播放、信号的重采样、信号加入噪声、信号的傅里叶变换和滤波等，以及GUI图形用户界面设计。 Ⅰ.语音信号的时域分析 语音信号是一种非平稳的时变信号，它携带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。语音信号分析的目的就在与方便有效的提取并表示语音信号所携带的信息。语音信号分析可以分为时域和变换域等处理方法，其中时域分析是最简单的方法。 Ⅱ.语音信号的频域分析 信号的傅立叶表示在信号的分析与处理中起着重要的作用。因为对于线性系统来说，可以很方便地确定其对正弦或复指数和的响应，所以傅立叶分析方法能完善地解决许多信号分析和处理问题。另外，傅立叶表示使信号的某些特性变得更明显，因此，它能更
MATLAB 在数字信号处理中的应用：连续信号的采样与重建 一、 设计目的和意义 随着通信技术的迅速发展以及计算机的广泛应用，利用数字系统处理模拟信号的情况变得更加普遍。数字电子计算机所处理和传送的都是不连续的数字信号，而实际中遇到的大都是连续变化的模拟量，现代应用中经常要求对模拟信号采样，将其转换为数字信号，然后对其进行计算处理，最好在重建为模拟信号。 采样在连续时间信号与离散时间信号之间其桥梁作用，是模拟信号数字化的第一个步骤，研究的重点是确定合适的采样频率，使得既要能够从采样信号（采样序列）中五失真地恢复原模拟信号，同时由要尽量降低采样频率，减少编码数据速率，有利于数据的存储、处理和传输。 本次设计中，通过使用用MATLAB 对信号f （t ）=A1sin(2πft)+A2sin(4πft)+A3sin(5πft)在300Hz 的频率点上进行采样，并进行仿真，进一步了解MA TLAB 在数字信号处理上的应用，更加深入的了解MA TLAB 的功能。 二、 设计原理 1、 时域抽样定理 令连续信号 xa(t)的傅立叶变换为Xa （j Ω），抽样脉冲序列p(t)傅立叶变换为P （j Ω），抽样后的信号x^(t)的傅立叶变换为X^(j Ω)若采用均匀抽样，抽样周期Ts ，抽样频率为Ωs= 2πfs ，有前面分析可知：抽样过程可以通过抽样脉冲序列p （t ）与连续信号xa （t ）相乘来完成，即满足：x^(t)p(t),又周期信号f （t ）傅立叶变换为： F[f(t)]=2[(]n s n F j n π δ∞ =-∞Ω-Ω∑ 故可以推得p(t)的傅立叶变换为： P （j Ω）=2[(]n s n P j n π δ∞ =-∞Ω-Ω∑ 其中： 根据卷积定理可知： X （j Ω）=12π Xa （j Ω）*P(j Ω) 得到抽样信号x （t ）的傅立叶变换为： X （j Ω）= [()]n n s n P X j n ∞=-∞Ω-Ω∑ 其表明：信号在时域被抽样后，他的频率X （j Ω）是连续信号频率X （j Ω）的形状以抽样频率Ωs 为间隔周期重复而得到，在重复过程中幅度被p （t ）的傅立叶级数Pn 加权。因为只是n 的函数，所以X （j Ω）在重复过程中不会使其形状发生变化。 假定信号x （t ）的频谱限制在-Ωm~+Ωm 的范围内，若以间隔Ts 对xa （t ）进行抽样信号X^（j Ω）是以Ωs 为周期重复。显然，若早抽样过程中Ωs<Ωm ，则 X^ (j Ω)将会发生频谱混叠的现象，只有在抽样的过程中满足Ωs>2Ωm 条件，X^(j Ω)才不会产生混频的混叠，在接收端完全可以有x^（t ）恢复原连续信号xa （t ），这就是低通信号的抽样定理的核心内容。
MATLAB中FFT的使用方法 说明：以下资源来源于《数字信号处理的MATLAB实现》万永革主编 一.调用方法 X=FFT(x)； X=FFT(x，N)； x=IFFT(X); x=IFFT(X,N) 用MATLAB进行谱分析时注意： （1）函数FFT返回值的数据结构具有对称性。 例： N=8; n=0:N-1; xn=[43267890]; Xk=fft(xn) → Xk= 39.0000-10.7782+6.2929i0-5.0000i 4.7782-7.7071i 5.0000 4.7782+7.7071i0+5.0000i-10.7782-6.2929i Xk与xn的维数相同，共有8个元素。Xk的第一个数对应于直流分量，即频率值为0。 （2）做FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。在IFFT时已经做了处理。要得到真实的振幅值的大小，只要将得到的变换后结果乘以2除以N即可。 二.FFT应用举例 例1：x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t)。采样频率fs=100Hz，分别绘制N=128、1024点幅频图。
clf; fs=100;N=128;%采样频率和数据点数 n=0:N-1;t=n/fs;%时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求得Fourier变换后的振幅 f=n*fs/N;%频率序列 subplot(2,2,1),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; subplot(2,2,2),plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅xlabel('频率/Hz'); ylabel('振幅');title('N=128');grid on; %对信号采样数据为1024点的处理 fs=100;N=1024;n=0:N-1;t=n/fs; x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t);%信号 y=fft(x,N);%对信号进行快速Fourier变换 mag=abs(y);%求取Fourier变换的振幅 f=n*fs/N; subplot(2,2,3),plot(f,mag);%绘出随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; subplot(2,2,4) plot(f(1:N/2),mag(1:N/2));%绘出Nyquist频率之前随频率变化的振幅 xlabel('频率/Hz'); ylabel('振幅');title('N=1024');grid on; 运行结果：
MATLAB 下的数字信号处理实现示例 附录一信号、系统和系统响应 1、理想采样信号序列 （1）首先产生信号x(n),0<=n<=50 n=0:50; A=444.128; a=50*sqrt(2.0)*pi; T=0.001; w0=50*sqrt(2.0)*pi; x=A*exp(-a*n*T).*sin(w0*n*T); close all subplot(3,1,1);stem(x); %定义序列的长度是50 %设置信号有关的参数 %采样率 %pi 是MATLAB 定义的π，信号乘可采用“.*”%清除已经绘制的x(n)图形 %绘制x(n)的图形 title(…理想采样信号序列?); （2）绘制信号x(n)的幅度谱和相位谱 k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n?*k); magX=abs(X); %绘制x(n)的幅度谱subplot(3,1,2);stem(magX);title(…理想采样信号序列的幅度谱?); angX=angle(X); %绘制x(n)的相位谱subplot(3,1,3);stem(angX) ; title (…理想采样信号序列的相位谱?) （3）改变参数为：A = 1,? = 0.4, & 0 = 2.0734, T = 1 n=0:50; A=1; a=0.4; T=1; w0=2.0734; x=A*exp(-a*n*T).*sin(w0*n*T); close all subplot(3,1,1);stem(x); title(…理想采样信号序列?); k=-25:25; W=(pi/12.5)*k; X=x*(exp(-j*pi/12.5)).^(n?*k); magX=abs(X); %定义序列的长度是50 %设置信号有关的参数 %采样率 %pi 是MATLAB 定义的π，信号乘可采用“.*” %清除已经绘制的x(n)图形 %绘制x(n)的图形 %绘制x(n)的幅度谱 subplot(3,1,2);stem(magX);title(…理想采样信号序列的幅度谱?); angX=angle(X); %绘制x(n)的相位谱
实验一 基于Matlab 的数字信号处理基本操作 一、 实验目的：学会运用MA TLAB 表示的常用离散时间信号；学会运用MA TLAB 实现离 散时间信号的基本运算。 二、 实验仪器：电脑一台，MATLAB6.5或更高级版本软件一套。 三、 实验内容： （一） 离散时间信号在MATLAB 中的表示 离散时间信号是指在离散时刻才有定义的信号，简称离散信号，或者序列。离散序列通常用)(n x 来表示，自变量必须是整数。 离散时间信号的波形绘制在MATLAB 中一般用stem 函数。stem 函数的基本用法和plot 函数一样，它绘制的波形图的每个样本点上有一个小圆圈，默认是空心的。如果要实心，需使用参数“fill ”、“filled ”，或者参数“.”。由于MATLAB 中矩阵元素的个数有限，所以MA TLAB 只能表示一定时间范围内有限长度的序列；而对于无限序列，也只能在一定时间范围内表示出来。类似于连续时间信号，离散时间信号也有一些典型的离散时间信号。 1. 单位取样序列 单位取样序列)(n δ，也称为单位冲激序列，定义为 ) 0() 0(0 1)(≠=?? ?=n n n δ 要注意，单位冲激序列不是单位冲激函数的简单离散抽样，它在n =0处是取确定的值1。在MATLAB 中，冲激序列可以通过编写以下的impDT .m 文件来实现，即 function y=impDT(n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 【实例1-1】 利用MATLAB 的impDT 函数绘出单位冲激序列的波形图。 解：MATLAB 源程序为 >>n=-3:3; >>x=impDT(n); >>stem(n,x,'fill'),xlabel('n'),grid on >>title('单位冲激序列') >>axis([-3 3 -0.1 1.1]) 程序运行结果如图1-1所示。 图1-1 单位冲激序列
工程设计论文 题目：基于MATLAB的语音信号采集与处理 姓名： 班级： 学号： 指导老师：
一.选题背景 1、实践意义： 语音信号是一种非平稳的时变信号，它携带着各种信息。在语音编码、语音合成、语音识别和语音增强等语音处理中无一例外需要提取语音中包含的各种信息。语音信号分析的目的就在于方便有效地提取并表示语音信号所携带的信息。所以理解并掌握语音信号的时域和频域特性是非常重要的。 通过语音相互传递信息是人类最重要的基本功能之一.语言是人类特有的功能.声音是人类常用工具,是相互传递信息的最重要的手段.虽然,人可以通过多种手段获得外界信息,但最重要,最精细的信息源只有语言,图像和文字三种.与用声音传递信息相比,显然用视觉和文字相互传递信息,其效果要差得多.这是因为语音中除包含实际发音容的话言信息外,还包括发音者是谁及喜怒哀乐等各种信息.所以,语音是人类最重要,最有效,最常用和最方便的交换信息的形式.另一方面,语言和语音与人的智力活动密切相关,与文化和社会的进步紧密相连,它具有最大的信息容量和最高的智能水平。 语音信号处理是研究用数字信号处理技术对语音信号进行处理的一门学科,处理的目的是用于得到某些参数以便高效传输或存储;或者是用于某种应用,如人工合成出语音,辨识出讲话者,识别出讲话容,进行语音增强等. 语音信号处理是一门新兴的学科,同时又是综合性的多学科领域,
是一门涉及面很广的交叉学科.虽然从事达一领域研究的人员主要来自信息处理及计算机等学科.但是它与语音学,语言学,声学,认知科学,生理学,心理学及数理统计等许多学科也有非常密切的联系. 语音信号处理是许多信息领域应用的核心技术之一,是目前发展最为迅速的信息科学研究领域中的一个.语音处理是目前极为活跃和热门的研究领域,其研究涉及一系列前沿科研课题,巳处于迅速发展之中;其研究成果具有重要的学术及应用价值. 数字信号处理是利用计算机或专用处理设备，以数值计算的方法对信号进行采集、抽样、变换、综合、估值与识别等加工处理，借以达到提取信息和便于应用的目的。它在语音、雷达、图像、系统控制、通信、航空航天、生物医学等众多领域都获得了极其广泛的应用。具有灵活、精确、抗干扰强、度快等优点。 数字滤波器, 是数字信号处理中及其重要的一部分。随着信息时代和数字技术的发展，受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。 FIR滤波器结构上主要是非递归结构，没有输出到输入的反馈，系统函数H (z)在处收敛，极点全部在z = 0处（因果系统），因而只能
短时能量分析matlab源程序： x=wavread('4.wav'); %计算N=50，帧移=50时的语音能量 s=fra(50,50,x);%对输入的语音信号进行分帧，其中帧长50，帧移50 s2=s.^2;%一帧内各种点的能量 energy=sum(s2,2);%求一帧能量 subplot(2,2,1); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=50'); axis([0,500,0,30]) %计算N=100，帧移=100时的语音能量 s=fra(100,100,x); s2=s.^2; energy=sum(s2,2); subplot(2,2,2); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=100'); axis([0,300,0,30]) %计算N=400，帧移=400时的语音能量 s=fra(400,400,x); s2=s.^2; energy=sum(s2,2); subplot(2,2,3); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=400'); axis([0,60,0,100]) %计算N=800，帧移=800时的语音能量 s=fra(800,800,x); s2=s.^2; energy=sum(s2,2); subplot(2,2,4); plot(energy) xlabel('帧数'); ylabel('短时能量E'); legend('N=800'); axis([0,30,0,200]) 分帧子函数： function f=fra(len,inc,x) %对读入语音分帧，len为帧长，inc为帧重叠样点数，x为输入语音数据 fh=fix(((size(x,1)-len)/inc)+1);%计算帧数 f=zeros(fh,len);%设一个零矩阵，行为帧数，列为帧长 i=1;n=1; while i<=fh %帧间循环 j=1; while j<=len %帧内循环 f(i,j)=x(n); j=j+1;n=n+1; end n=n-len+inc;%下一帧开始位置 i=i+1; end
大学本科毕业设计论文 基于MATLAB的有噪声语音信号处理
摘要 滤波器设计在数字信号处理中占有极其重要的地位,FIR数字滤波器和IＩR 滤波器是滤波器设计的重要组成部分。Matｌａｂ功能强大、简单易学、编程效率高,深受广大科技工作者的欢迎。特别是Matlab还具有信号分析工具箱,不需具备很强的编程能力,就可以很方便地进行信号分析、处理和设计。利用ＭATＬＡＢ信号处理工具箱可以快速有效地设计各种数字滤波器。课题基于MAＴLAB 有噪音语音信号处理的设计与实现，综合运用数字信号处理的理论知识对加噪声语音信号进行时域、频域分析和滤波。通过理论推导得出相应结论,再利用MＡTLAB作为编程工具进行计算机实现。在设计实现的过程中，使用窗函数法来设计FＩR数字滤波器，用巴特沃斯、切比雪夫和双线性变法设计IIＲ数字滤波器，并利用MATLＡB作为辅助工具完成设计中的计算与图形的绘制。通过对对所设计滤波器的仿真和频率特性分析，可知利用ＭＡＴLAB信号处理工具箱可以有效快捷地设计FIR和IIR数字滤波器，过程简单方便，结果的各项性能指标均达到指定要求。 关键词?数字滤波器；MATＬＡB；窗函数法；巴特沃斯; 切比雪夫; 双线性变换
Abｓtｒａｃt ?Filｔeｒdesiｇｎin ｄｉgｉtal sｉgnal proｃeｓsiｎｇplａysａn eｘtre ｍeｌｙiｍportａnt rｏle, FIR digital filtｅrs ａｎd IIR filter is an importａn ｔpaｒt of filter design.Matlab is ｐowｅrful,easy to ｌｅaｒn,pｒｏgｒａmming eｆfiｃieｎcｙ，which wａs wｅlcomｅd bｙthｅmajority oｆsc ｉentistｓ. Matlab alsｏhaｓ a particuｌar signaｌaｎalysis tｏolbｏx，it need notｈaｖe stronｇprｏgraｍminｇｓkiｌls can be easily signal anａlyｓiｓ, proｃesｓｉng aｎd design. Usｉng MATＬAB Signal Ｐroｃessing Toolbｏx can ｑｕickly anｄeｆficieｎtly design ａｖariｅtｙof digｉtaｌｆiltｅrs. ＭＡTLAB ｂaseｄon the noｉse issｕｅｓｐeech sｉｇnal processing deｓｉgn and ｉｍｐlｅmｅｎtation of digｉtal sｉgｎaｌprocessｉｎg inｔegrated use of ｔhe thｅoｒｅｔical knｏwlｅdge oｆtｈe ｓｐeecｈsｉｇｎal pｌus noiｓe, tｉme domaiｎ, frequencｙｄoｍaiｎanalysis ａnｄfiltering. Theｃorreｓponｄinｇresｕlts ｏbtａineｄthrougｈtheｏｒeticaｌｄerｉｖatｉon, ａnd thｅn use MATＬAB as a pｒogramｍing tooｌfor compuｔer iｍｐlementation．Iｍpleｍenｔｅd iｎthe desｉgn process，usinｇtｈｅｗindｏw function methoｄｔｏｄesiｇn FIR digiｔal fｉｌteｒs wｉｔh Butｔｅｒwortｈ, Chebyshev anｄbilinear Reｆｏｒm IIR dｉgｉtal filter dｅｓign ａｎd use ｏｆMATLAB ａs ａｓupplemenｔary tｏｏl to cｏｍｐｌete thｅcalculation and graphｉc dｅsiｇn Ｄrａwiｎg. Throughｔｈｅsimulａtion of ｔｈｅdesigned filter aｎd the fｒｅｑuenｃy anａlysis shoｗs thaｔuｓｉｎｇMatｌabＳｉｇnal Processing Toolbox can qｕickly and easily dｅsign digiｔal fｉlters FIR anｄIIR，the ｐroceｓｓis sｉｍple aｎd coｎvenｉent, the reｓuｌts of the perfｏrmaｎce iｎdｉcａtｏrs to ｍｅｅｔｔhe spｅｃｉfieｄrｅqｕirｅｍents. ? Kｅywｏrｄs: diｇiｔａl filｔer; ＭAＴLAB；Chebysｈev；Buttｅrｗｏrth；
实验一熟悉MATLAB环境 一、实验目的 (1)熟悉MATLAB的主要操作命令。 (2)学会简单的矩阵输入和数据读写。 (3)掌握简单的绘图命令。 (4)用MATLAB编程并学会创建函数。 (5)观察离散系统的频率响应。 二、实验内容 认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。 上机实验内容： (1)数组的加、减、乘、除和乘方运算。输入A=[1234]，B=[345 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出 A、B、C、D、E、F、G。 (2)用MATLAB实现以下序列。 a）x(n)=0.8n0≤n≤15 b）x(n)=e(0.2+3j)n0≤n≤15 c）x(n)=3cos(0.125πn+0.2π)+2sin(0.25πn+0.1π)0≤n≤15 (n)=x(n+16),绘出四个d）将c)中的x(n)扩展为以16为周期的函数x 16 周期。 (n)=x(n+10),绘出四个e）将c)中的x(n)扩展为以10为周期的函数x 10 周期。
(3)x(n)=[1,-1,3,5]，产生并绘出下列序列的样本。 a）x 1(n)=2x(n+2)-x(n-1)-2x(n) b）∑=-=5 1k 2) k n (nx (n) x (4)绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注。 a)x(t)=sin(2πt)0≤t≤10s b)x(t)=cos(100πt)sin(πt) 0≤t≤4s (5)编写函数stepshift(n0,n1,n2)实现u(n-n0),n1 数字信号处理（第二版） 绪论 1.4 MATLAB 在信号处理中的应用简介 MATLAB 是美国Mathworks 公司于1984年推出的一套高性能的数值计算和可视化软件，它集数值分析、矩阵运算、信号处理、系统仿真和图形显示于一体，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。 MATLAB 软件包括五大通用功能：数值计算功能（Numeric ），符号运算功能（Symbolic ）；数据可视化功能（Graphic ），数据图形文字统一处理功能（Notebook ）和建模仿真可视化功能（Simulink ）。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。目前，Mathworks 公司已推出30多个应用工具箱。MA TLAB 在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。 2.10 离散时间信号与系统的Matlab 表示 2.10.1 离散时间信号的表示和运算 1、基本序列的Matlab 表示 单位采样序列 在MA TLAB 中，单位采样序列可以通过编写以下的DTimpulse .m 文件来实现，即 function y=DTimpulse (n) y=(n==0); %当参数为0时冲激为1,否则为0 调用该函数时n 必须为整数或整数向量。 单位阶跃序列 在MA TLAB 中，单位阶跃序列可以通过编写DTu .m 文件来实现，即 function y=DTu (n) y=n>=0; %当参数为非负时输出1 调用该函数时n 必须为整数或整数向量。 矩形序列 用MA TLAB 表示矩形序列可根据公式()()()N R n u n u n N =--并利用DTu 函数生成，即 function y=DTR(n,N) y=DTu(n)-DTu(n-N); 调用该函数时n 必须为整数或整数向量，N 必须为整数。 实指数序列 用MA TLAB 表示实指数序列()(),n x n a u n n N a R =∈∈，即数字信号处理基本知识点Matlab实现
相关主题
文本预览