北理工数字信号处理实验报告.docx
《北理工数字信号处理实验报告.docx》由会员分享,可在线阅读,更多相关《北理工数字信号处理实验报告.docx(34页珍藏版)》请在冰豆网上搜索。
![北理工数字信号处理实验报告.docx](https://file1.bdocx.com/fileroot1/2023-2/3/45405341-7aae-439f-9620-1a02a0baa3d2/45405341-7aae-439f-9620-1a02a0baa3d21.gif)
北理工数字信号处理实验报告
数字信号处理实验报告
姓名:
仲易
班级:
05110901
学号:
20091087
实验1基2-FFT算法实现
一、实验目的
1.掌握基2-FFT的原理及具体实现方法。
2.编程实现基2-FFT算法。
3.深刻理解FFT算法的特点。
二、实验设备与环境
计算机,matlab软件环境。
三、实验基础理论
FFT是EDF的一种快速算法,能使DFT的计算大大简化,运算时间缩短。
FFT利用了WNnk的三个固有特性。
即对称性,周期性和可约性,将长序列的DFT分解为短序列的DFT,合并了DFT运算中的某些项,从而减少了DFT的运算量。
FFT算法基本上可分为两大类,即按时间抽取法和按频率抽取法。
在实现FFT算法时,要重点考虑两个问题,注意数据的读取和存储:
(1)输入输出的排序;
(2)蝶形运算的实现。
按时间抽取算法中输入反序输出顺序,按频率抽取算法中输入顺序输出反序;运算过程中的每一级都有N/2个蝶形运算构成,每一个蝶形运算单元中,两个节点变量运算后得到的结果为下一列相同位置的节点变量,而和其他节点变量无关,可以采用原位运算,节省存储单元。
另外,蝶形运算中的复系数WNnk可以存储为能及时查阅的系数表,这样可以借阅运算量,但是需要N/2哥复数存储器。
MATLAB中提供了用于计算FFT的函数fft,可将实验中所得到的结果与利用MATLAB中fft函数计算的结果相比较,以此验证结果的正确性。
4.实验内容及实验过程
1.编程实现序列长度为N=8的按时间抽取的基2-FFT算法。
给定一个8点序列,采用编写的程序计算其DFT,并与MATLAB中fft函数计算的结果相比较,以验证结果的正确性。
代码如下:
x=[12345678];
m=3;%X序列长度为2的3次幂
N=8;
nxd=bin2dec(fliplr(dec2bin([1:
N]-1,m)))+1;%求倒序列
y=x(nxd);
formm=1:
m%对x做m次基2分解
Nz=2^mm;u=1;
WN=exp(-i*2*pi/Nz);
forj=1:
Nz/2
fork=j:
Nz:
N
kp=k+Nz/2;
t=y(kp)*u;%蝶形运算乘项
y(kp)=y(k)-t;%加项
y(k)=y(k)+t;
end
u=u*WN;
end
end
y
y1=fft(x)%对比matlab的fft
y=
Columns1through3
36.0000-4.0000+9.6569i-4.0000+4.0000i
Columns4through6
-4.0000+1.6569i-4.0000-4.0000-1.6569i
Columns7through8
-4.0000-4.0000i-4.0000-9.6569i
y1=
Columns1through3
36.0000-4.0000+9.6569i-4.0000+4.0000i
Columns4through6
-4.0000+1.6569i-4.0000-4.0000-1.6569i
Columns7through8
-4.0000-4.0000i-4.0000-9.6569i
2、编程实现序列长度为N=8的按频率抽取的基2-FFT算法。
给定一个8点序列,采用编写的程序计算其DFT,并与MATLAB中fft函数计算的结果相比较,以验证结果的正确性。
代码如下:
L=2048;
M=256;
>>x=ones(1,L);
>>n=0:
M-1;
>>h=cos(0.2*pi*n);
>>tic
y1=conv(x,h);
toc
tic
X=fft(x,L+M-1);
H=fft(h,L+M-1);
Y=X.*H;
y2=ifft(Y);
toc
运行结果如下:
(分别为线形卷积和快速卷积)
Elapsedtimeis0.003352seconds.
Elapsedtimeis0.007167seconds.
L=4096;
M=256;
>>x=ones(1,L);
>>n=0:
M-1;
>>h=cos(0.2*pi*n);
>>tic
y1=conv(x,h);
toc
tic
X=fft(x,L+M-1);
H=fft(h,L+M-1);
Y=X.*H;
y2=ifft(Y);
toc
运行结果如下:
(分别为线形卷积和快速卷积)
Elapsedtimeis0.003666seconds.
Elapsedtimeis0.030134seconds.
3.将上述FFT程序推广大序列长度为N=2v的情况,要求利用原位运算。
代码如下:
L=2048;
M=256;
>>n=0:
M-1;
>>h=cos(0.2*pi*n);
l=128;
x=ones(1,l);
q=0;
tic
H=fft(h,l+M-1);
X=fft(x,l+M-1);
Y=X.*H;
y=ifft(Y);
yp=y;
fori=2:
16
q=q+l;
yq=[zeros(1,q),y];
yp=[yp,zeros(1,l)];
yp=yq+yp;
end
toc
运行结果如下:
Elapsedtimeis0.018875seconds.
4.
L=2048;
M=256;
n=0:
M-1;
h=cos(0.2*pi*n);
x=ones(1,L);
Lx=2048;
M=256;
M1=M-1;
N=512;
L=N-M1;
tic
h=[h,zeros(1,N-M)];
x=[zeros(1,M1),x,zeros(1,N-1)];
a=floor((Lx+M1-1)/(L))+1;
Y=zeros(1,N);
fork=0:
a-1
xk=x(k*L+1:
k*L+N);
b=fft(xk,N);
C=fft(h,N);
Z=b.*C;
Y(k+1,:
)=ifft(Z,N);
end
Y=Y(:
M:
N)';
Y=(Y(:
))'
toc
运行结果如下:
Elapsedtimeis0.035329seconds.
实验2利用DFT分析信号频谱
一、实验目的
1、加深对DFT原理的理解。
2、应用DFT分析信号的频谱。
3、深刻理解利用DFT分析信号频谱的原理,分析实现过程中出现的现象及解决方法。
二、实验设备与环境
计算机、MATLAB软件环境
三、实验基础理论
1、DFT与DTFT的关系
序列
的N点DFT
实际上就是
序列的DTFT在N个等间隔频率点
上的采样样本。
2、利用DFT求DTFT
方法一:
利用差值运算
其中
为内插函数
方法二:
在实际MATLAB计算中,采用增加数据的长度N,使得到的DFT谱线更加精细,用其包络近似计算DTFT。
如果数据量不足,可采用补零来增加数据长度。
3、利用DFT分析连续时间信号的频谱
1)确定时域采样间隔T,得到离散序列
2)确定截取长度M,得到M点离散序列
,这里
为窗函数
3)确定频域采样点数N,要求
4)利用FFT计算离散序列的N点DFT,得到
5)由
计算
采样点的近似值。
4、可能用到的MATLAB函数与代码
实验中DFT运算采用MATLAB中提供的函数fft来实现。
DTFT可以利用MATLAB矩阵运算的方法进行计算。
四、实验内容
1、
(1)
n=0:
3;
x=[2-111];
w=-pi:
0.01*pi:
pi;
X=x*exp(-j*n'*w);
subplot(211)
plot(w,abs(X));axistight
subplot(212)
plot(w,angle(X)/pi);axistight
(2)
在原基础上增加代码:
n=0:
3;
x=[2-111];
N=4;
k=0:
3;
Y=fft(x);
w=0:
0.01*pi:
2*pi;
X=x*exp(-j*n'*w);
m=0:
2*pi/4:
3/2*pi
subplot(211);
plot(w,abs(X));xlabel('\Omega/pi');title('DTFT-Magnitude');axistight
holdon
stem(m,abs(Y));
subplot(212);
plot(w,angle(X)/pi);xlabel('\Omega/pi');title('DTFT-Phase');
holdon
stem(m,angle(Y))
(3)
在原本基础上设置代码
subplot(211);
Z=fft(x,64);%n=64DFT
N=0:
63;
stem(abs(Z));
subplot(212);
stem(angle(Z))
2、考察序列
(1)
,用DFT估计
的频谱;将
补零加长到长度为100点序列用DFT估计
的频谱。
要求画出相应波形。
程序代码:
程序输出结果:
补零后:
程序代码:
程序输出结果:
(2)
时,用DFT估计
的频谱,并画出波形。
程序代码:
程序输出结果:
(3)根据实验结果,分析怎样提高频谱分辨率。
答:
减小截取长度或减小采样间隔的方法可以提高频谱分辨率。
3、已知信号
其中
=1Hz,
=2Hz,
=3Hz。
从
的表达式可以看出,它包含三个频率的正弦波,但是,从其时域波形来看,似乎是一个正弦信号,利用DFT做频谱分析,确定适合的参数,使得到的频谱的频率分辨率符合需要。
程序代码:
程序输出结果:
利用DFT近似分析连续时间信号
的频谱。
分析采用不同的采样间隔和截取长度进行计算的结果,并最终确定适合的参数。
程序代码:
程序输出结果:
通过改变步长可以达到改变采样间隔的结果,输出不同的频谱。
实验5、脉冲响应不变法设计IIR数字滤波器
一、实验目的
1、掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法
2、加深理解数字滤波器和模拟滤波器之间的技术指标转化。
3、掌握脉冲相应不变法设计IIR数字滤波器的优缺点及适用范围。
二、实验设备与环境
计算机、MATLAB软件环境
三、实验基础理论
1、基本原理:
从时域响应出发,使数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的单位冲激响应。
2、设计步骤:
a.确定数字滤波器性能指标。
b.将数字滤波器频率指标转换成相应的模拟滤波器的频率指标
c.设计模拟滤波器
实现方法:
脉冲响应不变法
[r,p,k]=residue(b,a)
[b,a]=residue(r,p,k)
d.模拟滤波器数字化
实现方法:
[bz,az]=impinvar(b,a,fs)
四、实验内容
1、设采样频率fs=4khz,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器,其3db截止频率为fc=1khz.
首先设计巴特沃斯滤波器:
N=3;
Wc=0.25*pi;
[b,a]=butter(N,Wc,'s')
b=
0000.4845
a=
1.00001.57081.23370.4845
因此,
然后脉冲响应不变法计算:
[bz,az]=impinvar(b,a)
bz=
1.1702-1.42120.6078-0.0945
az=
1.0000-1.62860.9286-0.1879
得到:
画出模拟滤波器和数字滤波器的图像进行对比:
[H,w]=freqs(b,a);
subplot(221);
plot(w/pi,abs(H));
gridon;
subplot(222);
plot(w/pi,20*log10((abs(H))/max(abs(H))));
w=[0:
500]*pi/500;
[H,w]=freqz(bz,az);
gridon;
subplot(223);
plot(w/pi,abs(H));
gridon;
subplot(224);
plot(w/pi,20*log10((abs(H))/max(abs(H))));
gridon;
2.设采样频率为fs=10khz,设计数字滤波器满足如下指标
通带截止频率:
1khz,通带波动:
1db,
阻带截止频率:
1.5khz,阻带衰减:
15db
分别采用巴特沃斯、切比雪夫I型、切比雪夫ii型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。
结合实验结果,分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求,分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。
(1)巴特沃斯滤波器为:
Wp=0.1*pi;
Ws=0.15*pi;
Rp=1;
As=15;
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws)))
OC1=Wp/((10^(Rp/10)-1)^(1/(2*N)))
[b,a]=butter(N,OC1,'s')
N=
6
OC1=
0.3516
b=
0000000.0019
a=
1.00001.35850.92270.39740.11410.02080.0019
得到:
[bz,az]=impinvar(b,a)
bz=
1.0e-003*
0.00000.01250.25730.51940.16370.00510
az=
1.0000-4.65179.1374-9.68295.8300-1.88890.2570
得到bz,az后,画出模拟滤波器和数字滤波器的图像进行对比:
[H,w]=freqs(b,a);
subplot(221);
plot(w/pi,abs(H));
gridon;
subplot(222);
plot(w/pi,20*log10((abs(H))/max(abs(H))));
w=[0:
500]*pi/500;
[H,w]=freqz(bz,az);
gridon;
subplot(223);
plot(w/pi,abs(H));
gridon;
subplot(224);
plot(w/pi,20*log10((abs(H))/max(abs(H))));
gridon;
[H,w]=freqs(b,a);
subplot(221);
plot(w/pi,abs(H));
gridon;
subplot(222);
plot(w/pi,20*log10((abs(H))/max(abs(H))));
w=[0:
500]*pi/500;
[H,w]=freqz(bz,az);
gridon;
subplot(223);
plot(w/pi,abs(H));
gridon;
subplot(224);
plot(w/pi,20*log10((abs(H))/max(abs(H))));
gridon;
(2)切比雪夫Ⅰ型
>>Wp=0.2*pi;
>>Ws=0.3*pi;
>>Rp=1;
>>As=15;
>>e=
>>e=(10^(Rp/10)-1)^(1/2)
e=0.5088
>>A=10^(As/20)
A=5.6234
>>N=ceil(acosh((A*A-1)^(1/2)/e)/acosh(Ws/Wp))N=4
>>[b,a]=cheby1(N,Rp,Wp,'s')
b=00000.0383a=1.00000.59870.57400.18420.0430
(3)切比雪夫Ⅱ型:
>>Wp=0.2*pi;
>>Ws=0.3*pi;
>>Rp=1;
>>As=15;
>>N=4;
>>[b,a]=cheby2(N,As,Ws,'s')
b=0.1778-0.00001.2637-0.00001.1225a=1.00002.36963.03221.99251.1225
(4)椭圆模拟原型滤波器:
>>Wp=0.2*pi;
>>Ws=0.3*pi;
>>Rp=1;
>>As=15;
>>e=0.5088;
>>A=5.6234;
>>k1=e/((A*A-1)^(1/2))
k1=0.0919
>>k=Wp/Ws
k=0.6667
>>[K,E]=ellipke(k)
K=2.0290
E=1.2612
>>[K1,E1]=ellipke((1-k1*k1)^(1/2))
K1=4.1217
E1=1.0077
>>[K2,E2]=ellipke(k1)
K2=1.6089
E2=1.5340
>>[K3,E3]=ellipke((1-k*k)^(1/2))
K3=2.1483
E3=1.2140
>>N=(K*K1)/(K2*K3)
N=2.4195
>>N2=2;
>>[b,a]=ellip(N2,Rp,As,Ws,'s')
b=0.177800.9379
a=1.00000.91201.0523
结果分析:
(1)当给定的设计指标相同时,选用椭圆滤波器所要求的结束N最低,切比雪夫滤波器
次之,巴特沃斯滤波器最高;
(2)从通带的相位相应来看,椭圆滤波器虽然提供了最优秀的幅度平方相应,但通带上
的相位相应非线性较大,而巴特沃斯滤波器在通带上具有相当的线性相位,切比雪
夫滤波器的相位特征介于两者之间。
实验七、窗函数法设计FIR数字滤波器
一、实验目的
掌握窗函数法设计FIR数字滤波器的原理及具体方法。
二、实验设备及环境
计算机、matlab软件及环境
三、实验基础理论
1.基本原理
窗函数设计法的基本思想为,首先选择一个适当的理想的滤波器Hd(ejw),然后用窗函数截取它的单位脉冲响应hd(n),得到线性相位和因果的FIR滤波器。
这种方法的重点是选择一个合适的窗函数和理想滤波器,使设计的滤波器的单位响应尽力逼向滤波器的单位响应。
2.设计步骤
(1)给定理想滤波器的频率响应Hd(ejw),再通带上具有单位增益和线性相位,在阻带上具有零响应。
(2)确定这个滤波器的单位脉冲响应
为了得到一个h(n)长度为N的因果相位FIR滤波器,我们令
(3)用窗函数截取hd(n)得到所设计FIR数字滤波器h(n)
h(n)=hd(n)w(n)
四、实验内容
1.设计一个数字低通滤波器,其技术指标如下
0.2
分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。
结合实验结果,分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。
2.设计一个数字低通滤波器,其技术指标如下
下阻带边缘:
下通带边缘:
0.35
上通带边缘:
0.65
上阻带边缘:
五、代码
确定长度:
wp=0.2*pi;
>>wst=0.3*pi;
>>width=wst-wp;
>>N=ceil(1.8*pi/width)+1
N=
19
确定hd(n)
>>n=0:
(N-1);
>>wc=(wp+wst)/2;
>>alpha=(N-1)/2;
>>hd=(wc/pi)*sinc((wc/pi)*(n-alpha))
>>wboxcar=boxcar(N)';
>>h=hd.*wboxcar
subplot(221);
>>stem(n,hd,'filled');
axistight;xlabel('n');ylabel('hd(n)');
>>[Hr,w1]=zerophase(h);
>>subplot(222);
plot(w1/pi,Hr);
>>axistight;xlabel('\omega^pi');ylabel('\omega');
>>subplot(223);
stem(n,h,'filled');
axistight;xlabel('n');ylabel('h(n)');
>>[H,w]=freqz(h,1);
>>subplot(224);
>>plot(w/pi,20*log10(abs(H)/max(abs(H))));
>>xlabel('\omega^pi');ylabel('dB')
>>gridon
由上图可清楚看出:
加矩形窗不满足指标
用海明窗:
wp=0.2*pi;
wst=0.3*pi;
width=wst-wp;
N=ceil(6.6*pi/width)+1
N=
67
n=0:
(N-1);
wc=(wp+wst)/2;
alpha=(N-1)/2;
hd=(wc/pi)*sinc((wc/pi)*(n-alpha));
whamming=hamming(N)';
h=hd.*whamming;
>>subplot(221);
stem(n,hd,'filled');
axistight;xlabel('n');ylabel('hd(n)');
[Hr,w1]=zerophase(h);
subplot(222);
plot(w1/pi,Hr);
axistight;xlabel('\omega^pi');ylabel('\omega');
subplot(223);
stem(n,h,'filled');
axistight;xlabel('n');ylabel('h(n)');
[H,w]=freqz(h,1);
subplot(224);
plot(w/pi,20*log10(abs(H)/max(abs(H))));
xlabel('\omega^pi');ylabel('dB')
gridon
由上图可清楚看出:
加海明窗满足指标
汉宁窗:
代码:
wp=0.2*pi;wst=0.3*pi;
tr_width=wst-wp;
N=ceil(6.2*pi/tr_width)+1;
n=0:
(N-1);wc=(wp+wst)/2;
alpha=(N-1)/2;
hd=(wc/pi)*sinc((wc/pi)*(n-alpha));
w_hanning=hanning(N)';
h=hd.*w_hanning;
subplot(221);stem(n,hd,'