实验指导书 36章.docx
《实验指导书 36章.docx》由会员分享,可在线阅读,更多相关《实验指导书 36章.docx(22页珍藏版)》请在冰点文库上搜索。
实验指导书36章
实验三离散傅里叶变换
一、实验目的
1.学习编制离散傅里叶变换程序.
2.理解DFT在数字信号处理中的核心地位和作用.
二、实验内容
1.编制计算离散傅里叶变换程序.
2.用编制程序处理时间抽样信号.
3.根据实序列离散傅里叶变换的对称性,初步判定程序的正确性.
三、实验说明
1.离散傅里叶变换公式如下
构造离散傅里叶正、反变换函数的MATLAB实现程序如下,其中dft(xn,N)为离散傅里叶正变换,idft(Xk,N)为离散傅里叶反变换:
---------------------------------------------------------------------------------------------------------------
function[Xk]=dft(xn,N)
n=[0:
1:
N-1];
k=n;
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk
function[xn]=idft(Xk,N)
n=[0:
1:
N-1];
k=n;
WN=exp(-j*2*pi/N);
nk=n'*k;
WNnk=WN.^(-nk);
xn=(Xk*WNnk)/N;
---------------------------------------------------------------------------------------------------------------
2.利用离散傅里叶变换函数求解序列傅里叶变换的MATLAB实现程序如下(N=8,
a=0.7)。
---------------------------------------------------------------------------------------------------------------
clear
N=8;
a=0.7;
n=[0:
7];
xn=a.^n;
Xk=dft(xn,N);
subplot(3,1,1)
stem(n,xn,'.k');axis([0,8,0,1.5])
subplot(3,1,2)
stem(n,abs(Xk),'.k');axis([0,8,0,5])
subplot(3,1,3)
stem(n,angle(Xk),'.k');axis([0,8,-1.5,1.5])
3.下图所示为
的幅度特性与相位特性。
其中,N为总的抽样点数;T为抽样时间间隔.
在这种条件下分析结果中的
以
点左右对称,说明该程序基本正确,可以进行其它分析.
4.在实验中同学们会发现,抽样信号分析得到的
的幅度与连续傅里叶变换分析该信号幅度不相同,这是值得大家深入讨论的问题.
3.14已知序列
=(3,11,7,0,-1,4,2),令
为
加入噪声干扰并移位后的序列:
其中,
为具有零均值和单位方差的高斯序列。
用MATLAB语言计算
和
之间的互相关。
3.14解:
程序及运行结果如下。
Clear
n=[1:
7];
xn=zeros(1,length(n));
xn2=xn;wn=xn;
xn(1:
7)=[3,11,7,0,-1,4,2];%signal
xn2(1:
7)=[4,2,3,11,7,0,-1];%delaysignal
wn(1:
7)=randn(1,7);%noise
yn=xn+wn;%signal+noise
xyn=conv(xn,yn);
subplot(2,2,1);
stem(n,xn,'k');title('x(n)');
subplot(2,2,2);
stem(n,xn2,'k');title('x(n-2)');
subplot(2,2,3);
stem(n,wn,'k');title('w(n)');
subplot(2,2,4);
stem(1:
(2*length(n)-1),xyn,'k');title('xy(n)');
四、实验报告要求
1.整理好经过运行并证明是正确的程序,并且加上详细的注释.
2.用连续傅里叶变换分析(被抽样的)连续信号,将其结果与抽样信号的离散傅里叶变换结果相比较,你能发现什么问题?
如何解释?
3.计算抽样序列的连续傅里叶变换,将其结果与抽样序列的离散傅里叶变换结果相比较,你又能发现什么问题?
如何解释?
五、实验结果分析
实验四快速傅里叶变换
一、实验目的
1.学习时间抽选奇偶分解FFT算法.
2.深入理解和掌握时间抽选奇偶分解FFT计算程序.
3.研究如何利用FFT程序分析确定性时间连续信号.
二、实验内容
1.用MATLAB编程比较DFT和FFT的运算时间
2.用MATLAB编程实现DFT和FFT的运算.
已知有限长序列
长度为N=4,
且:
1
2
–1
3
用FFT求X(k),再用IFFT求x(n).
3.用MATLAB程序分析FFT取不同长度时,序列x(n)的频谱变化情况.
三、实验说明
1.可以用以下MATLAB程序比较DFT和FFT的运算时间:
---------------------------------------------------------------------------------------------------------------
N=1024;
M=80;
x=[1:
M,zeros(1,N-M)];
t=cputime;
y1=fft(x,N);
Time_fft=cputime-t;
t1=cputime;
y2=dft(x,N);
Time_dft=cputime-t1;
t2=cputime;
---------------------------------------------------------------------------------------------------------------
Time_dft=
6.0290
Time_fft=
0.0100
由此可见FFT算法比直接计算DFT速度快得多。
2.利用快速傅里叶变换函数求解FFT和IFFT运算的MATLAB实现程序如下:
---------------------------------------------------------------------------------------------------------------
clear
xn=[1,2,-1,3];
X=fft(xn)
x=ifft(X)
---------------------------------------------------------------------------------------------------------------
X=[5.0000,2.0000+1.0000i,-5.0000,2.0000-1.0000i]
x=[1,2,-1,3]
3.设x(n)是长度为N=6的矩形序列,用MATLAB分析FFT取不同长度时x(n)的频谱变化。
N=8,32,64时x(n)的FFTMATLAB实现程序如下
---------------------------------------------------------------------------------------------------------------
x=[1,1,1,1,1,1];
N=8;y1=fft(x,N);
n=0:
N-1;subplot(3,1,1);stem(n,abs(y1),'.k');axis([0,9,0,6]);
N=32;y2=fft(x,N);
n=0:
N-1;subplot(3,1,2);stem(n,abs(y2),'.k');axis([0,40,0,6]);
N=64;y3=fft(x,N);
n=0:
N-1;subplot(3,1,3);stem(n,abs(y3),'.k');axis([0,80,0,6]);
---------------------------------------------------------------------------------------------------------------
x(n)的频谱如下图所示,N值取的越大就越接近序列真正的频谱。
4.8设
是长度为ML的长序列,其中
把
分成M段,记为
,
每段长度为L。
设
为L点冲激响应,则
显然,
是
点序列。
在这种方法中,需要保存中间卷积结果,在相加之前进行恰当的重叠,形成
。
(1)利用循环卷积,开发一个MATLAB函数实现重叠相加法;
(2)利用
(1)开发的函数采用基-2FFT,编写一个高速重叠相加分段卷积的MATLAB程序。
4.8解:
程序及运行结果如下。
(1)循环卷积子函数:
functionfn=circonvt(x1,x2,N);
%circonvt函数实现输入序列x1和x2的循环倦积,fn为输出序列
%N为循环卷积长度
if(length(x1)>N|length(x2>N)%判断输入信号的长度
error('N的长度必须大于输入数据的长度');
end
x1=[x1,zeros(1,N-length(x1))];
x2=[x2,zeros(1,N-length(x2))];
m=0:
N-1;
x=zeros(N,N);
forn=0:
N-1
x(:
n+1)=x2(mod((n-m),N)+1)';
end;
fn=x1*x;
主函数:
function[y]=ovrlpadd(x,h,L)
x=input('请输入x序列:
');
h=input('请输入y序列:
');
L=input('请输入段长L:
');
lenx=length(x);%x的长度
M=length(h);%h的长度
N1=L+M-1;%圆周卷积点数,即每一个输出序列Yi的长度
m=rem(lenx,L);%求余
ifm=0
x=[xzeros(1,L-m)];%末尾补零,使每段长度为N
K=floor(lenx/L)+1;%段数
else
x=x;
K=floor(lenx/L);
end
ytemp=zeros(1,N1-L);%N1-N为重叠部分,使其初始化为零
n1=1;n2=L;
fork=1;K
xk=x(n1:
n2);
Y(k,:
)=circonvt(xk,h,N1);
fori=1:
N1-L
Y(k,i)=Y(k,i)+ytemp(i);
ytemp(i)=Y(k,i+L);
end
y(n1:
n2+M-1)=Y(k,1:
N1);
n1=n1+L;n2=n2+L;
end
stem(y);
title('ovrlpadd');
请输入x序列:
[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
请输入y序列:
[1,0,0,-1]
请输入段长L:
4
ans=
1233333333333
333-14-15-16
(2)
function[y]=ovrlpaddfft(x,h,L)
x=input('请输入x序列:
');
h=input('请输入y序列:
');
L=input('请输入段长L:
');
lenx=length(x);%x的长度
M=length(h);%h的长度
N1=L+M-1;%满足循环卷积等于线性卷积的长度
N1=2^(ceil(log10(N1)/log10
(2)));
m=rem(lenx,L);
ifm=0
x=[xzeros(1,L-m)];%末尾补零,使每段长度为N
K=floor(lenx/L)+1;%段数
else
x=x;
K=floor(lenx/L);%段数
end
ytemp=zeros(1,N1-L);%N1-N为重叠部分,使其初始化为零
n1=1;
n2=L;
fork=1;K
z=x(n1:
n2)
xk=fft(z,N1);
Y(k,:
)=real(ifft(xk.*h));
fori=1:
N1-L
Y(k,i)=Y(k,i)+ytemp(i);
ytemp(i)=Y(k,i+L);
end
y(n1:
n2+N1-L)=Y(k,1:
N1);%输出结果
n1=n1+L;
n2=n2+L;
end
y=y(1:
(lenx+M-1));
stem(y);
title('ovrlpadd-2fft');
请输入x序列:
[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
请输入y序列:
[1,0,0,-1]
请输入段长L:
4
ans=
Columns1through12
1.00002.00003.00003.00003.00003.00003.00003.0000
3.00003.00003.00003.0000
Columns13through19
3.00003.00003.00003.0000-14.0000-15.0000-16.0000
四、实验报告要求
1.整理好经过运行并证明是正确的程序,并且加上详细的注释.
2.把N=1024点FFT分析所用的时间与直接计算DFT所用的时间相对比,从感性上理解”快”字.
3.深入讨论如何运用FFT来研究序列的频谱,为第八章的学习做好准备.
五、实验结果分析
实验五IIR滤波器设计
一.实验目的
1.学习模拟滤波器的设计方法.
2.学习模拟---数字变换滤波器设计方法.
3.掌握双线性变换数字滤波器设计法.
4.掌握用频带变换设计数字滤波器的具体方法.
二.实验内容
1.设计一个巴特沃斯模拟低通滤波器,满足以下性能指标:
通带的截止频率
,通带最大衰减
,阻带的截止频率
,阻带最大衰减
。
2.已知模拟滤波器的系统函数为
分别用冲激响应不变变焕法和双线性变换法将
转换为数字滤波器的系统函数
,并画出相应的
和
的频率响应曲线。
采样频率分别为
和
。
3.要求用双线性变换法设计一个数字巴特沃思低通滤波器,其特性曲线如下图所示。
在通带内
,允许幅度误差小于1dB,在阻带
时衰减应大于15dB。
通带幅度归一化,使其在
处为1。
4.用模拟频带变换法,由二阶巴特沃思函数设计截止频率为
抽样频率为
的数字高通滤波器。
三.实验说明
1.求解实验内容1的MATLAB实现程序如下:
---------------------------------------------------------------------------------------------------------------
clear;
closeall
fp=10000;fs=40000;Rp=3;As=35;
[N,fc]=buttord(fp,fs,Rp,As,'s')
[B,A]=butter(N,fc,'s');
[hf,f]=freqs(B,A,1024);
subplot(3,2,1);
plot(f,20*log10(abs(hf)/abs(hf
(1))))
grid;xlabel('f/Hz');ylabel('幅度(dB)');
axis([0,50000,-40,5])
line([0,50000],[-3,-3]);
---------------------------------------------------------------------------------------------------------------
程序运行结果:
N=3
fc=1.0441e+004
其频率特性曲线如下图所示:
频率单位:
2.求解实验内容2的MATLAB实现程序如下:
---------------------------------------------------------------------------------------------------------------
clear;
closeall
b=1000;a=[1,1000];
w=[0:
1000*2*pi];
[hf,w]=freqs(b,a,w);
subplot(2,3,1)
plot(w/2/pi,abs(hf));
grid;
xlabel('f/(Hz)');ylabel('幅度');
Fs0=[1000,500];
form=1:
2
Fs=Fs0(m);
[d,c]=impinvar(b,a,Fs);
wd=[0:
512]*pi/512;
hw1=freqz(d,c,wd);
subplot(2,3,2);
plot(wd/pi,abs(hw1)/abs(hw1
(1)));holdon;
end
grid;xlabel('f/(Hz)');
text(0.52,0.88,'T=0.002s');
text(0.12,0.54,'T=0.001s');
form=1:
2
Fs=Fs0(m);
[f,e]=bilinear(b,a,Fs);
wd=[0:
512]*pi/512;
hw2=freqz(f,e,wd);
subplot(2,3,3);
plot(wd/pi,abs(hw2)/abs(hw2
(1)));holdon;
end
grid;xlabel('f/(Hz)');
text(0.5,0.74,'T=0.002s');
text(0.12,0.34,'T=0.001s');
---------------------------------------------------------------------------------------------------------------
运行结果如下图所示:
f/Hz
(a)(b)(c)
模拟滤波器到数字滤波器的转换
由图(b)可见,对冲激响应不变法,采样频率越高(时间
越小),混叠越小;由图(c)可见,对双线性变换法,无频率混叠,但存在非线性失真。
3.求解实验内容3的MATLAB实现程序如下:
---------------------------------------------------------------------------------------------------------------
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
T=1;
Fs=1/T;
OmegaP=(2/T)*tan(wp/2);
OmegaS=(2/T)*tan(ws/2);
ep=sqrt(10^(Rp/10)-1);
Ripple=sqrt(1/(1+ep*ep));
Attn=1/(10^(As/20));
N=ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(OmegaP/OmegaS)));
OmegaC=OmegaP/((10^(Rp/10)-1)^(1/(2*N)));
[B,A]=butter(N,OmegaC,'s');
W=(0:
500)*pi/500;
[H]=freqs(B,A,W);
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
[b,a]=bilinear(B,A,T);
[h,w]=freqz(b,a,1000,'whole');
h=(h(1:
501))';
w=(w(1:
501))';
m=abs(H);
db=20*log10((m+eps)/max(m));
figure
(1);
subplot(2,2,1);plot(w/pi,mag);title('幅度')
ylabel(’模拟滤波器’);
axis([0,1,0,1.1])
set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]);
set(gca,’YTickmode’,’manual’,’YTick’,[0,Attn,Ripple,1]);grid
subplot(2,2,2);plot(w/pi,db1);title(’幅度(dB)’)
axis([0,1,-30,5])
set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]);
set(gca,’YTickmode’,’manual’,’YTick’,[-30,-15,-1,0]);grid
subplot(2,2,3);plot(w/pi,m);
xlabel(’频率单位:
pi’);ylabel(’数字滤波器’);
axis([0,1,0,1.1])
set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]);
set(gca,’YTickmode’,’manual’,’YTick’,[0,Attn,Ripple,1]);grid
subplot(2,2,4);plot(w/pi,db2);
xlabel(’频率单位:
pi’);
axis([0,1,-30,5])
set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]);
set(gca,’YTickmode’,’manual’,’YTick’,[-30,-15,-1,0]);grid
---------------------------------------------------------------------------------------------------------------
运行结果如下图所示。
N=6
OmegaC=0.7273
巴特沃思模拟滤波器及利用双线性变换法设计的数字滤波器
4.求解实验内容4的MATLAB实现程序如下:
---------------------------------------------------------------------------------------------------------------
N=2;
Fs=500;
fch=200;
wch=2*pi*fch/Fs;
[z,p,k]=buttap(N);
[b,a]=zp2tf(z,p,k);
[h,w]=freqs(b,a,512);
Omegach=2*