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

实验三 离散傅里叶变换

一、实验目的

1. 学习编制离散傅里叶变换程序.

2. 理解DFT在数字信号处理中的核心地位和作用. 二、实验内容

1. 编制计算离散傅里叶变换程序. 2. 用编制程序处理时间抽样信号.

3. 根据实序列离散傅里叶变换的对称性,初步判定程序的正确性. 三、实验说明

1.离散傅里叶变换公式如下

N?1X(k)? ??n?0N?1n?0x(n)e?j2?Nnk?x(n)cos??2???2??nk??j?x(n)sin?nk??N??N?n?0N?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.下图所示为X?k?的幅度特性与相位特性。

其中,N为总的抽样点数;T为抽样时间间隔. 在这种条件下分析结果中的X?k?,以其它分析.

4.在实验中同学们会发现,抽样信号分析得到的X?k?的幅度与连续傅里叶变换分析该信号幅度不相同,这是值得大家深入讨论的问题.

N2点左右对称,说明该程序基本正确,可以进行

3.14 已知序列x(n)=(3,11,7,0,-1,4,2),令y(n)为x(n)加入噪声干扰并移位后的序列：

y(n)?x(n-2)??(n)

其中,?(n)为具有零均值和单位方差的高斯序列。用MATLAB语言计算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]; Tlay signal 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的运算. 已知有限长序列x(n)长度为N ＝4, 且:

1 n?0 x(n)? 2 n?1 –1 n?2

3 n?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)的FFT MATLAB实现程序如下 ---------------------------------------------------------------------------------------------------------------

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 设x?n?是长度为ML的长序列，其中M?1,L?1,把x?n?分成M段，记为xm?n?，

m?1,2,?,M,每段长度为L。

?x?n?xm?n????0

mM?n??m?1?M?1其他M?1

x?n???x?n?

mM?1m?0设h(n)为L点冲激响应，则

y?n??x?n??h?n??M?1?x?n??h?n???mym?n?

m?0m?0ym?n??xm?n??h?n?

显然，ym?n?是2L?1点序列。在这种方法中，需要保存中间卷积结果，在相加之前进行恰当的重叠，形成y?n?。

(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的长度必须大于输入数据的长度'); end

x1=[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-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); %求余 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 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);

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 =

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=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);

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 Columns13through19

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.设计一个巴特沃斯模拟低通滤波器，满足以下性能指标：通带的截止频率

?p?10000rad/s，通带最大衰减Ap?3dB，阻带的截止频率?s?40000rad/s，阻

带最大衰减As?35dB。

2. 已知模拟滤波器的系统函数为

Ha(s)?1000s?1000

分别用冲激响应不变变焕法和双线性变换法将Ha(s)转换为数字滤波器的系统函数

H(z)，并画出相应的Ha(s)和H(z)的频率响应曲线。采样频率分别为1000Hz和500Hz。

3. 要求用双线性变换法设计一个数字巴特沃思低通滤波器，其特性曲线如下图所示。在通带内?p?0.2?，允许幅度误差小于1dB，在阻带?d?0.3?时衰减应大于15dB。通带幅度归一化，使其在??0处为1。

4. 用模拟频带变换法，由二阶巴特沃思函数设计截止频率为200Hz,抽样频率为

500Hz的数字高通滤波器。

三.实验说明

1.求解实验内容1的MATLAB实现程序如下：

---------------------------------------------------------------------------------------------------------------

clear;

close all

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

其频率特性曲线如下图所示:

频率 单位：10000rad/s

2.求解实验内容2的MATLAB实现程序如下：

---------------------------------------------------------------------------------------------------------------

clear; close all

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]; 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=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）可见，对冲激响应不变法，采样频率越高（时间T越小），混叠越小；由图(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)));求解截止频率?c [B,A]=butter(N,OmegaC,'s');巴特沃斯滤波器设计 W=(0:500)*pi/500; W取值范围和抽样间隔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'); 采用单位圆上1000个点 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('幅度') 4个图中的图一，命名幅度 ylabel(’模拟滤波器’)；y坐标轴名 axis([0,1,0,1.1]) 坐标轴设置x.y轴的范围

set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]); x轴取样点坐标

set(gca,’YTickmode’,’manual’,’YTick’,[0,Attn,Ripple,1]);grid y轴取样点坐标 subplot(2,2,2);plot(w/pi,db1);title(’幅度（dB）’) 四个字图中的图二，命名幅度 axis([0,1,-30,5]) 坐标轴设置x.y轴的范围

set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]); x轴取样点坐标

set(gca,’YTickmode’,’manual’,’YTick’,[-30,-15,-1,0]);grid y轴取样点坐标

subplot(2,2,3);plot(w/pi,m); 4个子图中的图三

xlabel(’频率单位：pi’);ylabel(’数字滤波器’)；坐标轴命名 axis([0,1,0,1.1]) 坐标轴设置x.y轴的范围

set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]); x轴取样点坐标

set(gca,’YTickmode’,’manual’,’YTick’,[0,Attn,Ripple,1]);grid y轴取样点坐标

subplot(2,2,4);plot(w/pi,db2);四个字图中的图四 xlabel(’频率单位:pi’);横坐标轴名 axis([0,1,-30,5]) 坐标轴设置x.y轴的范围

set(gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]); x轴取样点坐标

set(gca,’YTickmode’,’manual’,’YTick’,[-30,-15,-1,0]);grid y轴取样点坐标

--------------------------------------------------------------------------------------------------------------- 运行结果如下图所示。

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*Fs*tan(wch/2); [Bs,As]=lp2hp(b,a,Omegach); [Bz,Az]=bilinear(Bs,As,Fs) [H,W]=freqz(Bz,Az,512); m=abs(H);

db2=20*log10((m+eps)/max(m)); figure(1);

subplot(2,2,1);plot(w/pi,db1) ylable(’幅度（dB）’); axis([0,4,-40,0]);grid subplot(2,2,2);plot(W/pi,db2); ylable(‘幅度（dB）’);

axis([0,1,-100,0]);grid

---------------------------------------------------------------------------------------------------------------

运行结果如下：

Bz = 0.0675 -0.1349 0.0675 Az = 1.0000 1.1430 0.4128 其频带变换曲线如下图所示:

(a)模拟低通滤波器 (b)数字高通滤波器

低通到高通的频带变换曲线

四.实验报告要求

1.给出详细的滤波器设计说明书.

2. 整理好经过运行并证明是正确的程序,并且加上详细的注释. 3.验证所设计滤波器的幅度特性和相位特性.

4.用所设计的滤波器对不同的正弦信号进行滤波,以说明其特性. 5.对数字滤波器模拟截止频率与抽样时间的关系进行讨论. 五.实验结果分析

本文来源：https://www.bwwdw.com/article/zh3p.html