数字信号处理实验二Word文档格式.docx
《数字信号处理实验二Word文档格式.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验二Word文档格式.docx(30页珍藏版)》请在冰豆网上搜索。
![数字信号处理实验二Word文档格式.docx](https://file1.bdocx.com/fileroot1/2023-1/30/1c5e2f12-e6cb-4c04-bb7a-a42612a623f9/1c5e2f12-e6cb-4c04-bb7a-a42612a623f91.gif)
四、实验原理
参考《数字信号处理》教材
《数字信号处理的MATLAB实现》万永革编著
五、主要实验仪器及材料
微型计算机、Matlab。
六、实验步骤
1.设计一简单正弦信号,通过改变采样率观察采样前后的信号变化。
例如:
假设有一振幅为1,频率为10Hz,相位为0.3的模拟信号,即
,用0.01s的采样间隔(采样频率为100Hz)来表示原始信号(注意:
实际上模拟信号不能用离散值表示,此处为了在计算机上表示,用采样率非常高的离散信号表示模拟信号)。
分别以5Hz,10Hz(每秒采样10次,即采样间隔为0.1s),20Hz,40Hz,80Hz,200Hz对原始信号进行采样,画出采样前后的信号,并画出其频谱图,对比前后的变化,验证采样定理。
(1)可以用t=0:
1/fs:
9/f;
取9个周期,通过改变采样率,自动改变采样点数。
(2)也可以通过设置dt1(采样间隔),已知采样点数n1,t1=n1*dt1,
如图所示,采样率为40Hz时的原始信号,采样过程和采样后的信号时域图和频谱图,可见,当采样率大于原始信号频率的两倍时,采样前后信号频率基本不发生变化,信号不失真。
图2-1采样的过程
图2-2原始信号和采样后信号的频谱
2.语音信号的采集。
利用windows下的录音机或其他软件,录制一段自己的话音,时间控制在1秒左右。
然后在MATLAB软件平台下,利用函数wavread对语音信号进行采样,记住采样频率和采样点数。
通过wavread函数的使用,要求理解采样频率、采样位数等概念。
wavread函数调用格式:
y=wavread(file),读取file所规定的wav文件,返回采样值放在向量y中。
[y,fs,nbits]=wavread(file),采样值放在向量y中,fs表示采样频率(Hz),nbits表示采样位数。
y=wavread(file,N),读取前N点的采样值放在向量y中。
y=wavread(file,[N1,N2]),读取从N1点到N2点的采样值放在向量y中。
3.对真实语音信号的采样、重构。
(1)降采样:
利用windows下的录音机录制的音频采样率是固定的fs(=22050),可以选择以下函数实现对语音信号的降采样。
y=x((1:
N:
length(x)));
%对原始信号每隔N个点取一位,即采样率变为原来的1/N
y=resample(yn,L,M);
%采样率变为原来的L/M倍
y=downsample(yn,N);
%%采样率变为原来的1/N倍
改变采样率为原来的1/2倍,1/4倍,1/20倍,1/50倍,1/100倍等,分别画出降采样前后的信号波形和频谱图,分析采样前后信号的变化。
图2-3采样率为原来的1/2时的信号波形频谱图
图2-4采样率为原来的1/10时信号波形频谱图
在MATLAB中,函数sound可以对声音进行回放。
其调用格式:
sound(x,fs,bits);
通过调用此函数,感觉采样前后声音的变化。
(2)重构原信号:
降采样后,信号的采样率和采样点数同时变化。
如采样率变为原来的1/2,即对原始信号每隔一个点采样。
如果要恢复原始信号,即信号长度和采样率须变为原来同样大小。
所以,必须对降采样后信号进行插值重构。
具体过程参见附件1中例子。
对采样后的真实语音信号进行插值重构,滤波,恢复原始信号。
画出插值前后信号的波形以及频谱图,并将重构后信号与原始信号进行比较。
如,对采样率降为原来1/5的降采样后信号插值重构,结果如下图所示。
图2-5采样率为原来的1/5时的波形频谱图
图2-6插值后的信号波形频谱图
图2-7低通滤波器的频率响应图
图2-8滤波后的波形频谱图
调用sound函数感受插值后的声音,发现会有高频的噪声。
经过低通滤波器之后,高频噪声被滤出。
但是,因为之前的原始信号经过了降采样,所以插值后的效果一定不如原始声音。
也可用用wavwrite函数将经过处理的语音信号保存下来,调用格式为如wavwrite(y,fs,bits,'
sound.wav'
),其中,y为所要保存的语音信号,fs为其采样率,bits为采样位数,'
是保存的语音信号的文件名。
4.语音信号的频谱分析
要求:
首先画出语音信号的时域波形;
然后对语音信号进行频谱分析,在MATLAB中,可以利用函数fft对信号进行快速付立叶变换,得到信号的频谱特性;
从而加深对频谱特性的理解。
[x,fs,bits]=wavread('
e:
\sound1.wav'
[10245120]);
%读取[10245120]段
%语音信号数据;
如果不指定所读取读取语音信号数据的长度,则读取整段语音信号
X=fft(x,4096);
%fftshift(fft(x))
N=length(x);
k=0:
N-1;
subplot(221);
plot(k,x);
title('
原始信号波形'
);
subplot(222);
plot((-N/2:
N/2-1)/N*fs,abs(X));
原始信号频谱'
5.对原始语音信号加噪声:
设计频率已知的噪声信号或者用自然噪声信号加在原始语音信号上,构建带噪声信号。
如对原始信号加上频率为5000的正弦波噪声信号,程序如下:
sound1.wav'
%读取原始语音信号
X=abs(fftshift(fft(x)));
%原始信号频谱
c1=0.01*sin(2*pi*5000*k/fs);
%构建频率为5000Hz的正弦波噪声信号
yn=x+c1'
;
%构建带噪声信号
Yn=abs(fftshift(fft(yn)));
%求带噪声信号频谱
subplot(321);
plot(x);
subplot(322);
N/2-1)/length(k)*fs,X);
subplot(323);
plot(yn);
带噪声信号波形'
subplot(324);
N/2-1)/length(k)*fs,Yn);
带噪声信号频谱'
运行程序,结果如图2-9所示。
从频谱图可以看出,原始信号频谱段集中在0~5000Hz之间,主要频率集中在0Hz附近的低频部分。
加上噪声信号后,在5000Hz处有个幅值非常大的的高频成份,即以上所加的高频正弦波噪声信号。
图2-9加噪声前后信号的波形频谱图
6.设计数字滤波器和画出频率响应
根据分析所得的原始信号的频谱和噪声信号频谱特点,给出有关滤波器的性能指标。
首先用窗函数法或者最优化法设计高通,低通,带通,带阻滤波器,在MATLAB中,可以利用函数fir1,firls设计FIR滤波器;
然后在用双线性变换法或脉冲响应法设计上面几种滤波器,在MATLAB中,可以利用函数butte、cheby1和ellip设计IIR滤波器;
最后,利用MATLAB中的函数freqz画出各滤波器的频率响应。
具体方法参加附件3种滤波器设计的步骤和实例。
如,根据以上语音信号的特点给出有关IIR滤波器的性能指标:
1)低通滤波器性能指标,fp=4500,fc=6500,Rs=100,Rp=1。
(fp:
通带截至频率;
fc:
阻带截至频率;
Rs:
通带波纹;
Rp:
阻带波纹)
2)带阻滤波器性能指标,fp1=4800Hz,fp2=5200Hz,fc1=4600Hz,fc2=5400Hz,Rs=30dB,Rp=1dB。
([fp1fp2]:
[fc1fc2]:
通带截至频率)
程序如下:
%椭圆带阻滤波器
figure,
wp=[48005200]*2/fs;
ws=[46005400]*2/fs;
%通带和阻带边界频率
Rp=1;
Rs=30;
Nn=128;
%通带波纹和阻带衰减以及绘制频率特性的数据点数
[NN,Wn]=ellipord(wp,ws,Rp,Rs);
%求取数字滤波器的最小阶数和归一化截至频率
[b,a]=ellip(NN,Rp,Rs,Wn,'
stop'
%设计滤波器
freqz(b,a,512,fs);
title('
椭圆带阻滤波器'
图2-10椭圆带阻滤波器的频率响应
%双线性变换法设计的椭圆低通滤波器:
fp=4500;
fc=6500;
Rs=100;
fs=22050;
wc=2*fc/fs;
wp=2*fp/fs;
[n,wn]=ellipord(wp,wc,Ap,As);
[b,a]=ellip(n,Rp,Rs,wn);
若设计FIR滤波器,要阻止5000Hz的高频信号通过,即设计如下滤波器特性,
%最优化设计法设计带阻滤波器
n=100;
%滤波器阶数
f=[04600*2/fs4800*2/fs5400*2/fs5600*2/fs1];
%频率向量
a=[110011];
%振幅向量
b=firls(n,f,a);
%采用firls设计滤波器
[h,w1]=freqz(b);
%计算其频率响应
plot(w1/pi,abs(h));
%绘制滤波器的幅频响应
最优化法设计的带阻滤波器'
图2-11最优化法设计的带阻滤波器的幅频响应
7.用滤波器对信号进行滤波
比较各种滤波器的性能,然后用性能好的各滤波器分别对采集的信号进行滤波。
比较滤波前后语音信号的波形及频谱,要求在一个窗口同时画出滤波前后的波形及频谱。
在MATLAB中,FIR滤波器利用函数fftfilt对信号进行滤波,IIR滤波器利用函数filter对信号进行滤波。
xx=filter(b,a,y);
%对信号y滤波
XX=fftshift(fft(xx));
%求经过滤波后信号的傅立叶变换
图2-12对带噪声信号滤波
8、图形用户界面的设计
为使编制的程序操作方便,有兴趣的同学,可以利用Matlab进行图形用户界面的设计。
在所设计的系统界面上可以选择滤波器的类型,输入滤波器的参数,显示滤波器的频率响应,选择信号等。
七、实验报告要求
1.按照实验指导书的格式简述实验目的,基本要求,实验内容以及实验原理。
2.详细写明每一个实验的具体步骤,如何实现,在具体编程中遇到的问题以及如何解决。
3.按照实验要求编程,给出详细程序实现代码,以及程序结果(图,表格等)。
4.在每一步中,分析各种不同的情况所得出的结果,并进行对比,得出结论。
5.在本次实验报告的最后,总结本次实验的主要内容,以及所掌握的内容。
附件1:
信号的插值重构
以对一个简单正弦信号的重构为例,以10倍原信号长度进行插值。
clear;
closeall;
f=10;
fs=100;
k=1:
100;
x=sin(2*pi*f*k/fs);
n=10;
y=zeros(1,n*length(x));
y(1:
n:
length(y))=x;
附图2-1信号插值前后的波形
附图2-2信号插值前后的频谱图
由上图可看出,插值后信号最大频率同样变为原来的10倍,高频部分进行了拓展(由于在原始信号中插入0的缘故),为原始频谱的周期延拓(类似采样定理)。
故,只需根据以上插值后信号的频率特性,设计FIR低通滤波器进行滤波,即可得到(滤波器具体的设计方法见附录)。
低频滤波器的通带为0~50Hz,由于1对应最大频率500Hz,故,50Hz对应的转换频率为50/500=0.1,设计通带截止频率为wp为0.09,阻带截止频率为0.11。
N为滤波器阶数,理想状况下,阶数越高,滤波器效果越好。
(工程中涉及到的滤波器阶数问题需另做考虑)
N=100;
B=firls(N,[00.090.111],[1100]);
figure,plot((-N/2:
N/2)/(N/2),abs(fftshift(fft(B)))),title(‘FIR低通滤波器’);
xx=filter(B,1,y);
%利用低通滤波器滤波
subplot(211),plot(xx);
滤波后信号'
subplot(212),plot((-n*50:
n*50-1)/(n*100)*n*fs,abs(fftshift(fft(xx))));
滤波后信号频谱'
附图2-3FIR低通滤波器的频率响应
附图2-4插值后的信号经低通滤波的信号波形和频谱图
附件2:
信号频谱
有限长序列的离散傅立叶变换(DFT)定义为
逆变换为:
Matlab频谱分析相关函数:
(1)一维快速傅立叶变换FFT函数fft,调用格式:
y=fft(x)
y=fft(x,n)
函数说明:
fft函数用于计算矢量或矩阵的离散傅立叶变换,可通过
来实现,式中
。
Y=fft(x)利用FFT算法计算矢量x的离散傅立叶变换,当x为矩阵时,y为矩阵x的每一列的FFT。
当x的长度为2的幂次方时,则fft采用基-2FFT算法,否则采用稍慢的混合基算法。
Y=fft(x,n)采用n点FFT。
当x长度小于n时,fft函数在x尾部补零.。
当x长度大于n时,函数将序列截断.
(2)一维逆快速傅立叶变换(IFFT)ifft,调用格式:
Y=ifft(x)
Y=ifft(x,n)
ifft函数用于计算矢量或矩阵的逆傅立叶变换,即
式中,
(3)fft函数最通常的应用是计算信号的频谱。
考虑由一个50Hz和120Hz正弦信号构成的信号,受零均值随机信号的干扰,数据采样率为1000Hz。
通过fft函数来分析其信号频率成分。
相关程序如下:
t=0:
0.001:
0.6;
x=sin(2*pi*50*t)+sin(2*pi*120*t);
y=x+1.5*randn(1,length(t));
Y=fft(y,512);
附件3:
滤波器设计
一、IIR数字滤波器设计
1、经典设计法:
主要有两种方法:
脉冲响应不变法、双线性变换法
经典设计IIR滤波器的步骤:
A.根据给定的性能指标和方法,首先对设计性能指标中的频率指标进行转换,转换后的频率指标作为模拟滤波器原型设计指标;
B.估计模拟滤波器最小阶数和边界频率,可利用Matlab工具函数buttord,cheb1ord等;
如,设计buttord滤波器的调用格式:
[n,Wn]=buttord(Wp,Ws,Rp,Rs)
[n,Wn]=buttord(Wp,Ws,Rp,Rs,’s’)
buttord可在给定滤波器性能的情况下,选择模拟或数字滤波器的最小阶数,其中,Wp和Ws分别是通带和阻带的截止频率,取值范围为0<
Wp(或Ws)<
1,其值为1时表示0。
fs为采样频率,Rp和Rs分别是通带和阻带区的波纹系数。
C.设计模拟低通滤波器原型,可利用buttap,cheblap等;
D.由模拟低通原型经频率变换得到模拟滤波器(低通,高通,带通,带阻),可利用MATLAB工具函数lp2lp,lp2hp,lp2bs,lp2bp等;
如低通到低通模拟滤波器变换函数为[bt,at]=lp2lp(b,a,Wo),lp2lp函数将截止频率为1rad/s(归一化频率)的模拟低通滤波器原型变换为截止频率为Wo的低通滤波器。
E.将模拟滤波器离散化得到IIR数字滤波器,可利用MATLAB工具函数bilinear,impinvar等。
例1.用脉冲响应不变法设计Butterworth低通滤波器,要求通带频率为
,通带波纹小于1dB,阻带在
内,幅度衰减大于15dB,采样周期T=0.01s。
假设一个信号
,其中
试将原信号与经过该滤波器的输出信号进行比较。
注意,该题给出的通带边界频率和阻带边界频率均为数字频率,因此设计时首先要将其转换为模拟频率。
%Samp1
wp=0.2*pi;
ws=0.3*pi;
Rs=15;
%数字滤波器截止频率、通带波纹和阻带衰减
T=0.01;
%采样间隔
Wp=wp/T;
Ws=ws/T;
%得到模拟滤波器的频率—采用脉冲响应不变法的频率转换形式
[N,Wn]=buttord(Wp,Ws,Rp,Rs,'
s'
%计算模拟滤波器的最小阶数
[z,p,k]=buttap(N);
%设计低通原型数字滤波器
[Bap,Aap]=zp2tf(z,p,k);
%零点极点增益形式转换为传递函数形式
[b,a]=lp2lp(Bap,Aap,Wn);
%低通滤波器频率转换
[bz,az]=impinvar(b,a,1/T);
%脉冲响应不变法设计数字滤波器传递函数
figure
(1)
[H,f]=freqz(bz,az,Nn,1/T);
%输出幅频响应和相频响应
subplot(2,1,1),plot(f,20*log10(abs(H)));
xlabel('
频率/Hz'
ylabel('
振幅/dB'
gridon;
subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))
相位/^o'
figure
(2)
f1=5;
f2=30;
%输入信号含有的频率
%数据点数
n=0:
t=n*T;
%时间序列
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t);
%输入信号
subplot(2,1,1),plot(t,x),title('
输入信号'
)
y=filtfilt(bz,az,x);
%对信号进行滤波
subplot(2,1,2),plot(t,y),title('
输出信号'
),xlabel('
时间/s'
程序运行结果如下图所示。
该例所要求的通带频率为
,而该滤波器的采样间隔为0.01s,采样频率为100Hz,即
对应于100Hz,则
对应于10Hz,即该滤波器的通带范围为0~10Hz。
观看附图2-5,其通带范围最大衰减小于1dB,与该分析一致。
题意要求阻带在
内,幅度衰减大于15dB。
其中,
对应于15Hz,幅度衰减大于15dB。
完全符合滤波器设计的要求。
附图2-5例1所设计滤波器的幅频响应和相频响应
附图2-6例1所设计滤波器的输入输出信号
例2.用双线性变换法设计一个椭圆低通滤波器,其性能指标同例1。
%Samp2
%数字滤波器截止频率通带波纹和阻带衰减
Fs=100;
Ts=1/Fs;
%采样频率
Wp=2/Ts*tan(wp/2.);
Ws=2/Ts*tan(ws/2.);
%按频率转换公式进行转换——双线性变换法
[N,Wn]=ellipord(Wp,Ws,Rp,Rs,'
[z,p,k]=ellipap(N,Rp,Rs);
%设计模拟原型滤波器
%低通转换为低通滤波器的频率转换
[bz,az]=bilinear(b,a,Fs);
%运用双线性变换法得到数字滤波器传递函数
[H,f]=freqz(bz,az,Nn,Fs);
%求出频率特性
程序运行如下图,在10Hz以前,衰减小于1dB,在15Hz以后,衰减大于15dB,即性能指标完全满足滤波器设计的要求。
附图2-7例2所设计滤波器的幅频响应和相频响应
2、IIR滤波器完全设计函数:
以上介绍了IIR滤波器的基本方法和步骤,并给出一些例子说明如何利用MATLAB编程分步实现这些步骤。
由这些步骤可知,我们必须多次调用MATLAB信号处理工具箱中的基本工具函数。
实际上,MATLAB信号处理工具箱还提供了IIR滤波器设计的完全工具函数,用户只需要调用这些工具函数即可一次性地完成设计,而不需要调用那些基本工具函数分步实现。
数字IIR滤波器的完全设计函数有
[b,a]=butter(n,Wn,’ftype’)
[b,a]=cheby1(n,Rp,Wn,’ftype’)
[b,a]=cheby2(n,Rs,Wn,’ftype’)
[b,a]=ellip(n,Rp,Rs,Wn,’ftype’)
在上面的调用函数中,n代表滤波器阶数,Wn代表滤波器的截止频率,取值为0~1