实验指导书 36章.docx

1、实验指导书 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

2、.nk; Xk=xn*WNnkfunction 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)。-clearN=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,

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 解：程序及运行结果如下。 Clearn=1:7;xn=zeros(

4、1,length(n);xn2=xn;wn=xn;xn(1:7)=3,11,7,0,-1,4,2; %signalxn2(1:7)=4,2,3,11,7,0,-1; %delay signalwn(1:7)=randn(1,7); %noiseyn=xn+wn; %signal+noisexyn=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

5、(1:(2*length(n)-1),xyn,k);title(xy(n); 四、实验报告要求1.整理好经过运行并证明是正确的程序,并且加上详细的注释.2.用连续傅里叶变换分析(被抽样的)连续信号,将其结果与抽样信号的离散傅里叶变换结果相比较,你能发现什么问题?如何解释?3.计算抽样序列的连续傅里叶变换,将其结果与抽样序列的离散傅里叶变换结果相比较,你又能发现什么问题?如何解释? 五、实验结果分析实验四 快速傅里叶变换一、 实验目的1. 学习时间抽选奇偶分解FFT算法.2. 深入理解和掌握时间抽选奇偶分解FFT计算程序.3. 研究如何利用FFT程序分析确定性时间连续信号.二、实验内容1. 用M

6、ATLAB编程比较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;-Ti

7、me_dft = 6.0290Time_fft =0.0100由此可见FFT算法比直接计算DFT速度快得多。2.利用快速傅里叶变换函数求解FFT和 IFFT运算的MATLAB实现程序如下：-clearxn=1,2,-1,3;X=fft(xn)x=ifft(X)-X = 5.0000,2.0000 + 1.0000i,-5.0000,2.0000 - 1.0000ix = 1,2,-1,33.设x(n)是长度为N = 6的矩形序列，用MATLAB分析FFT取不同长度时x(n)的频谱变化。N=8，32，64时x(n)的FFT MATLAB实现程序如下-x=1,1,1,1,1,1;N=8;y1=ff

8、t(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点冲激响应，则显然，是点序列。在

9、这种方法中，需要保存中间卷积结果，在相加之前进行恰当的重叠，形成。(1) 利用循环卷积，开发一个MATLAB函数实现重叠相加法；(2) 利用（1）开发的函数采用基-2FFT，编写一个高速重叠相加分段卷积的MATLAB程序。4.8 解：程序及运行结果如下。(1)循环卷积子函数： function fn=circonvt(x1，x2，N)； %circonvt函数实现输入序列x1和x2的循环倦积，fn为输出序列 %N为循环卷积长度 if(length(x1)N|length(x2N)%判断输入信号的长度 error(N的长度必须大于输入数据的长度); end x1=x1,zeros(1,N-len

10、gth(x1); x2=x2,zeros(1,N-length(x2); m=0:N-1; x=zeros(N,N); for n=0:N-1 x(:,n+1)=x2(mod(n-m),N)+1); end; fn=x1*x;主函数: functiony=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); %求余 if m=0 x=x

11、 zeros(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; for k=1;K xk=x(n1:n2); Y(k,:)=circonvt(xk,h,N1); for i=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); tit

12、le(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 = 1 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 -14 -15 -16(2) functiony=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(

13、log10(N1)/log10(2); m=rem(lenx,L); if m=0 x=x zeros(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; for k=1;K z=x(n1:n2) xk=fft(z,N1); Y(k,:)=real(ifft(xk. *h); for i=1:N1-L Y(k,i)=Y(k,i)+ytemp(i); ytemp(i)=Y(k,i+L);

14、 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.0000 2.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 Columns13thr

15、ough19 3.0000 3.0000 3.0000 3.0000 -14.0000 -15.0000 -16.0000四、实验报告要求1.整理好经过运行并证明是正确的程序,并且加上详细的注释.2.把N=1024点FFT分析所用的时间与直接计算DFT所用的时间相对比,从感性上理解”快”字.3.深入讨论如何运用FFT来研究序列的频谱,为第八章的学习做好准备.五、实验结果分析实验五 IIR滤波器设计一. 实验目的1. 学习模拟滤波器的设计方法.2. 学习模拟-数字变换滤波器设计方法.3. 掌握双线性变换数字滤波器设计法.4. 掌握用频带变换设计数字滤波器的具体方法.二. 实验内容1.设计一个巴特

16、沃斯模拟低通滤波器，满足以下性能指标：通带的截止频率，通带最大衰减，阻带的截止频率，阻带最大衰减。2. 已知模拟滤波器的系统函数为 分别用冲激响应不变变焕法和双线性变换法将转换为数字滤波器的系统函数，并画出相应的和的频率响应曲线。采样频率分别为和。3. 要求用双线性变换法设计一个数字巴特沃思低通滤波器，其特性曲线如下图所示。在通带内，允许幅度误差小于1dB，在阻带时衰减应大于15dB。通带幅度归一化，使其在处为1。4. 用模拟频带变换法，由二阶巴特沃思函数设计截止频率为,抽样频率为的数字高通滤波器。三.实验说明1.求解实验内容1的MATLAB实现程序如下：-clear;close allfp=

17、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 = 3fc = 1.0441e+004其频率特性曲线如下图所示: 频率 单位： 2.求解实验内容2的MATLAB实现程序如下：-clear;close al

18、lb=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;for m=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);hold on;end grid;xlabel(f/(Hz); text(0.52,0.88,T

19、=0.002s); text(0.12,0.54,T=0.001s);for m=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);hold on;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）可见，对冲激响应不

20、变法，采样频率越高（时间越小），混叠越小；由图(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=Om

21、egaP/(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,

22、0,1.1)set(gca,XTickMode,manual,XTick,0,0.2,0.3,1);set(gca,YTickmode,manual,YTick,0,Attn,Ripple,1);gridsubplot(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);gridsubplot(2,2,3);plot(w/pi,m);xlabel(频率单位：pi);yla

23、bel(数字滤波器)；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);gridsubplot(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 = 6OmegaC = 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*