实验指导书 3-6章(1)

实验三 离散傅里叶变换一、实验目的1. 学习编制离散傅里叶变换程序.2. 理解DFT 在数字信号处理中的核心地位和作用. 二、实验内容1. 编制计算离散傅里叶变换程序.2. 用编制程序处理时间抽样信号.3. 根据实序列离散傅里叶变换的对称性,初步判定程序的正确性. 三、实验说明1.离散傅里叶变换公式如下21011()()22()cos ()sin N jnkNn N N n n X k x n ex n nk j x n nk N N πππ--=--===⎛⎫⎛⎫=-⎪ ⎪⎝⎭⎝⎭∑∑∑构造离散傅里叶正、反变换函数的MA TLAB 实现程序如下，其中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*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.利用离散傅里叶变换函数求解序列傅里叶变换的MA TLAB 实现程序如下(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,1,3)stem(n,angle(Xk),'.k');axis([0,8,-1.5,1.5]) 3.下图所示为()k X 的幅度特性与相位特性。

其中,N 为总的抽样点数;T 为抽样时间间隔. 在这种条件下分析结果中的()k X ,以2N 点左右对称,说明该程序基本正确,可以进行其它分析.4.在实验中同学们会发现,抽样信号分析得到的()k X 的幅度与连续傅里叶变换分析该信号幅度不相同,这是值得大家深入讨论的问题.3.14 已知序列()x n =(3,11,7,0,-1,4,2),令()y n 为()x n 加入噪声干扰并移位后的序列：()(-2)()y n x n n ω=+其中,()n ω为具有零均值和单位方差的高斯序列。

用MA TLAB 语言计算()x n 和()y n 之间的互相关。

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]; %delay signal wn(1:7)=randn(1,7); %noise yn=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(1:(2*length(n)-1),xyn,'k');title('xy(n)');四、实验报告要求1.整理好经过运行并证明是正确的程序,并且加上详细的注释.2.用连续傅里叶变换分析(被抽样的)连续信号,将其结果与抽样信号的离散傅里叶变换结果相比较,你能发现什么问题?如何解释?3.计算抽样序列的连续傅里叶变换,将其结果与抽样序列的离散傅里叶变换结果相比较,你又能发现什么问题?如何解释?五、实验结果分析实验四 快速傅里叶变换一、 实验目的1. 学习时间抽选奇偶分解FFT 算法.2. 深入理解和掌握时间抽选奇偶分解FFT 计算程序.3. 研究如何利用FFT 程序分析确定性时间连续信号. 二、实验内容1. 用MA TLAB 编程比较DFT 和FFT 的运算时间2. 用MA TLAB 编程实现DFT 和FFT 的运算. 已知有限长序列)(n x 长度为N ＝4, 且:1 0=n =)(n x2 1=n –1 2=n3 3=n 用FFT 求X （k ），再用IFFT 求x(n).3. 用MA TLAB 程序分析FFT 取不同长度时，序列x(n)的频谱变化情况. 三、实验说明1.可以用以下MA TLAB 程序比较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.0290Time_fft =0.0100由此可见FFT 算法比直接计算DFT 速度快得多。

2.利用快速傅里叶变换函数求解FFT 和 IFFT 运算的MA TLAB 实现程序如下： ---------------------------------------------------------------------------------------------------------------clearxn=[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的矩形序列，用MA TLAB分析FFT取不同长度时x(n)的频谱变化。

N=8，32，64时x(n)的FFT MA TLAB实现程序如下--------------------------------------------------------------------------------------------------------------- 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 设()n x 是长度为ML 的长序列，其中1,1M L ,把()n x 分成M 段，记为()n x m ，,,,2,1M m =每段长度为L 。

()()⎩⎨⎧=0n x n x m()其他11-+≤≤M m n mM()()∑-==1M m mn x n x设()h n 为L 点冲激响应，则()()()()()()∑∑-=-==*=*=110M m m M m mn y n h n x n h n x n y()()()n h n x n y m m *=显然，()n y m 是21L -点序列。

在这种方法中，需要保存中间卷积结果，在相加之前进行恰当的重叠，形成()n y 。

(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(x2>N)%判断输入信号的长度error('N的长度必须大于输入数据的长度');endx1=[x1,zeros(1,N-length(x1))];x2=[x2,zeros(1,N-length(x2))];m=0:N-1;x=zeros(N,N);for n=0:N-1x(:,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); %求余if m=0x=[x zeros(1,L-m)]; %末尾补零，使每段长度为NK=floor(lenx/L)+1; %段数elsex=x;K=floor(lenx/L);endytemp=zeros(1,N1-L); %N1-N为重叠部分，使其初始化为零n1=1;n2=L;for k=1;Kxk=x(n1:n2);Y(k,:)=circonvt(xk,h,N1);for i=1:N1-LY(k,i)=Y(k,i)+ytemp(i);ytemp(i)=Y(k,i+L);endy(n1:n2+M-1)=Y(k,1:N1);n1=n1+L;n2=n2+L;endstem(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: 4ans =1 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 -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);if m=0x=[x zeros(1,L-m)]; %末尾补零，使每段长度为NK=floor(lenx/L)+1; %段数elsex=x;K=floor(lenx/L); %段数endytemp=zeros(1,N1-L); %N1-N为重叠部分，使其初始化为零n1=1;n2=L;for k=1;Kz=x(n1:n2)xk=fft(z,N1);Y(k,:)=real(ifft(xk. *h));for i=1:N1-LY(k,i)=Y(k,i)+ytemp(i);ytemp(i)=Y(k,i+L);endy(n1:n2+N1-L)=Y(k,1:N1); %输出结果n1=n1+L;n2=n2+L;endy=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: 4ans =Columns1through121.00002.00003.0000 3.0000 3.0000 3.0000 3.0000 3.00003.0000 3.0000 3.0000 3.0000Columns13through193.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.设计一个巴特沃斯模拟低通滤波器，满足以下性能指标：通带的截止频率10000/p r a d s Ω=，通带最大衰减dB A p 3=，阻带的截止频率s rad s /40000=Ω，阻带最大衰减dB A s 35=。