数字信号处理实验fft进行谱分析，可以应用在什么方面

数字信号处理综合设计 

一、实验目的

1．学会MATLAB的使用，掌握MATLAB的程序设计方法；

2．掌握在Windows环境下语音信号采集的方法；

3．掌握数字信号处理的基本概念、基本理论和基本方法；

4．掌握MATLAB设计FIR和IIR数字滤波器的方法；

5．学会用MATLAB对信号进行分析和处理。

二、实验原理

参考《数字信号处理》教材。

三、主要实验仪器及材料

微型计算机、Matlab65教学版、TC编程环境。

四、实验内容

1．语音信号的采集

要求利用windows下的录音机（开始—程序—附件—娱乐—录音机，文件—属性—立即转换—8000KHz，8位，单声道）或其他，录制一段自己的话音，时间控制在1秒左右。然后在MATLAB下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解采样频率、采样位数等概念。

wavread函数调用格式：

y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。

[y,fs,nbits]=wavread(file)，采样值放在向量y中，fs表示采样频率（Hz），nbits表示采样位数。

y=wavread(file,N)，读取前N点的采样值放在向量y中。

y=wavread(file,[N1,N2])，读取从N1点到N2点的采样值放在向量y中。

2．语音信号的频谱分析

要求首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在MATLAB中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。

3．设计数字滤波器和画出频率响应

根据语音信号的特点给出有关滤波器的性能指标：1）低通滤波器性能指标，fp=1000Hz，fc=1200 Hz， As=100dB，Ap=1dB；2）高通滤波器性能指标，fc=2800 Hz，fp=3000 Hz As=100dB，Ap=1dB；3）带通滤波器性能指标，fp1=1200 Hz，fp2=3000 Hz，fc1=1000 Hz，fc2=3200 Hz，As=100dB，Ap=1dB。要求学生首先用窗函数法设计上面要求的三种滤波器，在MATLAB中，可以利用函数fir1设计FIR滤波器，然后在用双线性变换法设计上面要求的三种滤波器；之后再利用函数butter和cheby1设计上面要求的三种IIR滤波器。最后，利用MATLAB中的函数freqz画出各滤波器的频率响应。

4．用滤波器对信号进行滤波

比较FIR和IIR两种滤波器的性能，然后用性能好的各滤波器分别对采集的信号进行滤波，在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。

5．比较滤波前后语音信号的波形及频谱

要求在一个窗口同时画出滤波前后的波形及频谱。

6．回放语音信号

在MATLAB中，函数sound可以对声音进行回放。其调用格式：

sound(x,fs,bits)；

可以感觉滤波前后的声音有变化。

五、实验思考

1．双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？ 

2．能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？

六、实验报告要求 

简述实验原理及目的。

2按照实验步骤及要求，比较各种情况下的滤波性能。

3总结实验所得主要结论。

4简要回答思考题。

不是没有图，而是后3个画的只有一个点，查看这几个变量的值发现它们都是一个标量而不是向量或矩阵。一直往前推查，发现可能是你的s变量重复赋值了一次，变成了一个变量。

这是编程中很容易排查的小问题。你应该学会设置断点，并查看此时workspace中的变量，计算机是个很诚实的东西。

我试了一下，把第14行和17行的s换个名字应该能解决。（你的代码好像贴了两遍。。。敢不敢不要这么粗心。。。读的'speech'文件我还是随便找了个文件凑上去的）

s1=10^(-rs/20);

fpts=[wp ws];

mag=[0 1];

dev=[p s1];

答：数字频率度量不相同，但他们所对应的模拟频率相同。由w=ΩTs公式得，采样间隔变化时模拟频率对应的数字频率会有相应的变化，故其度量会有所变化。而且采样频率的大小直接关系到能否将能否将原始信号恢复出来。

直接有函数调用的滤波器分析和实现 看matlab的help文件，然后写个m文件就可以了

fftfilt 基于FFT重叠加法的数据滤波

filter 递归（IIR）或非递归（FIR）滤波器的数据滤波

firter2 二维数字滤波

filtfilt 零相位数字滤波

filtic 函数filter初始条件确定

freqs 模拟滤波器频率响应

freqspace 频率响应的频率空间设置

freqz 数字滤波器频率响应

impz 数字滤波器的脉冲响应

latcfilt 格型梯形滤波器实现

实验一 系统响应及系统稳定性

一、实验目的

（1）掌握 求系统响应的方法。

（2）掌握时域离散系统的时域特性。

（3）分析、观察及检验系统的稳定性。

二、实验内容

（1）给定一个低通滤波器的差分方程为

y(n)=005x(n)+005x(n-1)+09y(n-1)

输入信号x1(n)=R8(n)

x2(n)=u(n)

(a) 分别求出系统对x1(n)=R8(n) 和x2(n)=u(n)的响应序列，并画出其波形。

(b) 求出系统的单位冲响应，画出其波形。

实验程序：

A=[1,-09];B=[005,005]; %%系统差分方程系数向量 B 和 A

x1n=[1 1 1 1 1 1 1 1 zeros(1,50)]; %产生信号 x1(n)=R8(n)

x2n=ones(1,128); %产生信号 x2(n)=u(n)

y1n=filter(B,A,x1n); %求系统对 x1(n)的响应 y1(n)

n=0:length(y1n)-1;

subplot(2,2,1);

stem(n,y1n,"");

title("(a) 系统对 R_8(n)的响应y_1(n)");

xlabel("n");

ylabel("y_1(n)");

y2n=filter(B,A,x2n); %求系统对 x2(n)的响应 y2(n)

n=0:length(y2n)-1;

subplot(2,2,2);

stem(n,y2n,"");

title("(b) 系统对 u(n)的响应y_2(n)");

xlabel("n");

ylabel("y_2(n)");

hn=impz(B,A,58); %求系统单位脉冲响应 h(n)

n=0:length(hn)-1;

subplot(2,2,3);

y=hn;

stem(n,hn,"");

title("(c) 系统单位脉冲响应h(n)");

xlabel("n");

ylabel("h(n)");

运行结果图：

（2）给定系统的单位脉冲响应为

h1(n)=R10(n)

h2(n)= δ(n)+25δ(n-1)+25δ(n-2)+δ(n-3)

用线性卷积法分别求系统h1(n)和h2(n)对x1(n)=R8(n)的输出响应，波形。

实验程序：

x1n=[1 1 1 1 1 1 1 1 ]; %产生信号 x1(n)=R8(n)

h1n=ones(1,10);

h2n=[1 25 25 1 ];

y21n=conv(h1n,x1n);

y22n=conv(h2n,x1n);

figure(2)

n=0:length(h1n)-1;

subplot(2,2,1);

stem(n,h1n);

title("(d) 系统单位脉冲响应h1n");

xlabel("n");

ylabel("h_1(n)");

n=0:length(y21n)-1;

subplot(2,2,2);

stem(n,y21n);

title("(e) h_1(n)与 R_8(n)的卷积y_{21}n");

并画出

xlabel("n");

ylabel("y_{21}(n)");

n=0:length(h2n)-1;

subplot(2,2,3);

stem(n,h2n);

title("(f) 系统单位脉冲响应h_2n");

xlabel("n");

ylabel("h_2(n)");

n=0:length(y22n)-1;

subplot(2,2,4)

stem(n,y22n);

title("(g) h_2(n)与 R_8(n)的卷积y_{22}n");

xlabel("n");

ylabel("y_{22}(n)");

运行结果图：

（3）给定一谐振器的差分方程为

y(n)=18237y(n-1)-09801y(n-2)+b0x(n)-b0x(n-2)

令b0=1/10049，谐振器的谐振频率为 04rad。

(a) 用实验方法检查系统是否稳定。输入信号为u(n) 时，画出系统输出波形。

(b) 给定输入信号为

x(n)= sin(0014n )+ sin(04n )

求出系统的输出响应，并画出其波形。

实验程序：

un=ones(1,256); %产生信号 u(n)

n=0:255;

xsin=sin(0014n)+sin(04n); %产生正弦信号

A=[1,-18237,09801];

B=[1/10049,0,-1/10049]; %系统差分方程系数向量 B和 A

y31n=filter(B,A,un); %谐振器对 u(n)的响应 y31(n)

y32n=filter(B,A,xsin); %谐振器对 u(n)的响应 y31(n)

figure(3)

n=0:length(y31n)-1;

subplot(2,1,1);

stem(n,y31n,"");

title("(h) 谐振器对 u(n) 的响应y_{31}n");

xlabel("n");

ylabel("y_{31}(n)");

n=0:length(y32n)-1;

subplot(2,1,2);

stem(n,y32n,"");

title("(i) 谐振器对正弦信号的响应y_{32}n");

xlabel("n");

ylabel("y_{32}(n)");

运行结果图：

实验二 时域采样与频域采样

一、实验目的

时域采样理论与频域采样理论是数字信号处理中的重要理论。要求掌握模拟信号

采样前后频谱的变化，以及如何选择采样频率才能使采样后的信号不丢失信息；要求掌握频率域采样会引起时域周期化的概念，以及频率域采样定理及其对频域采样点数选择的指导作用。

二、实验内容

（1）时域采样理论的验证

给定模拟信号,x a (t)=Ae-at sin(Ω0t)u(t)

式中 A=444128，α =50 2 π，Ω =50 2 πrad/s

Tp=64/1000; %观察时间 Tp=64 微秒

Fs=1000;

T=1/Fs;

M=TpFs;

n=0:M-1;

t=nT;

A=444128;

alph=pi502^05;

omega=pi502^05;

xat=Aexp(-alpht)sin(omegat);%给定模拟信号

Xk=Tfft(xat,M); %M 点 FFT[xat)]

subplot(3,2,1); stem(n,xat,"");

xlabel("n");

ylabel("x_1(n)");

title("(a) Fs=1000Hz");

k=0:M-1;fk=k/Tp;

subplot(3,2,2);

plot(fk,abs(Xk));

title("(a) TFT[x_a(nT)],F_s=1000Hz");

xlabel("\omega/hz");

ylabel("(H_1(e^j^w))");

axis([0,Fs,0,12max(abs(Xk))]);

Fs=300;T=1/Fs;

M=TpFs;

n=0:M-1;

t=nT;

A=444128;

alph=pi502^05;

omega=pi502^05;

xat=Aexp(-alpht)sin(omegat);

Xk=Tfft(xat,M); %M 点 FFT[xat)]

subplot(3,2,3);

stem(n,xat,"");

xlabel("n");

ylabel("x_2(n)");

title("(b) Fs=300Hz");

k=0:M-1;fk=k/Tp;

subplot(3,2,4);

plot(fk,abs(Xk));

title("(a) TFT[x_a(nT)],Fs=300Hz");

xlabel("\omega/hz");

ylabel("(H_2(e^j^w))");

axis([0,Fs,0,12max(abs(Xk))]);

A=444128;

alph=pi502^05;

omega=pi502^05;

xat=Aexp(-alpht)sin(omegat);

Xk=Tfft(xat,M); %M 点 FFT[xat)]

subplot(3,2,5);

stem(n,xat,"");

xlabel("n");

ylabel("x_3(n)");

title("(c) Fs=200Hz");

k=0:M-1;fk=k/Tp;

subplot(3,2,6);

plot(fk,abs(Xk));

title("(a) TFT[x_a(nT)],Fs=200Hz");

xlabel("\omega/hz");

ylabel("(H_3(e^j^w))");

axis([0,Fs,0,12max(abs(Xk))])

（2）频域采样理论的验证

clc;clear;close all;

M=27;N=32;n=0:M;

xa=0:(M/2);

xb=ceil(M/2)-1:-1:0;

xn=[xa,xb]; %产生 M 长三角波序列 x(n)

Xk=fft(xn,1024); %1024点FFT[x(n)], 用于近似序列 x(n)的 TF X32k=fft(xn,32) ;%32 点 FFT[x(n)]

x32n=ifft(X32k); %32点 IFFT[X32(k)]得到 x32(n)

X16k=X32k(1:2:N); %隔点抽取 X32k 得到 X16(k)

x16n=ifft(X16k,N/2); %16点 IFFT[X16(k)]得到 x16(n)

subplot(3,2,2);

stem(n,xn,"");

title("(b) 三角波序列 x(n)");

xlabel("n");

ylabel("x(n)");

axis([0,32,0,20]);

k=0:1023;wk=2k/1024;

subplot(3,2,1);

plot(wk,abs(Xk));

title("(a)FT[x(n)]");

xlabel("\omega/\pi");

ylabel("|X(e^j^\omega)|");

axis([0,1,0,200]) ;

k=0:N/2-1;

subplot(3,2,3);

stem(k,abs(X16k),"");

title("(c) 16 点频域采样");

xlabel("k");

ylabel("|X_1_6(k)|");

axis([0,8,0,200])

n1=0:N/2-1;

subplot(3,2,4);

stem(n1,x16n,"");

title("(d) 16点 IDFT[X_1_6(k)]");

xlabel("n");

ylabel("x_1_6(n)");

axis([0,32,0,20]);

k=0:N-1;

subplot(3,2,5);

stem(k,abs(X32k),"");

title("(e) 32 点频域采样");

xlabel("k");

ylabel("|X_3_2(k)|");

axis([0,16,0,200]);

n1=0:N-1;

subplot(3,2,6);

stem(n1,x32n,"");

box on

title("(f) 32 点IDFT[X_3_2(k)]");

xlabel("n");

ylabel("x_3_2(n)");

axis([0,32,0,20]);

运行结果图：

实验三 用 FFT 对信号作频谱分析

一．实验目的

学习用FFT 对连续信号和时域离散信号进行谱分析的方法，了解可能出现

的分析误差及其原因，以便正确应用 FFT。

二．实验内容

（1）对以下序列进行谱分析

⎧n +10≤n ≤3⎪x 1(n )=R 4(n )；x 2(2)=⎨8-n 4≤n ≤7；

⎪0其它n ⎩

0≤n ≤3⎧4-n ⎪x 3(n ) =⎨n -34≤n ≤7

⎪0其它n ⎩

选择 FFT 的变换区间 N为 8 和16 两种情况进行频谱分析。分别打印其幅频特性曲线。 并进行对比、分析和讨论。

实验程序：

x1n=[ones(1,4)]; %产生序列向量x1(n)=R4(n)

M=8;xa=1:(M/2);

xb=(M/2):-1:1;

x2n=[xa,xb]; %产生长度为8 的三角波

x3n=[xb,xa];

X1k8=fft(x1n,8); %计算x1n 的 8点DFT

X1k16=fft(x1n,16); %计算x1n 的16点 DFT

X2k8=fft(x2n,8); %计算x2n 的 8点DFT

X2k16=fft(x2n,16); %计算x2n 的16点 DFT

X3k8=fft(x3n,8); %计算x3n 的 8点DFT

X3k16=fft(x3n,16); %计算x3n 的16点 DFT

%以下绘制幅频特性曲线

subplot(1,2,1);

stem(X1k8,""); %绘制8点 DFT 的幅频特性图

title("(1a) 8点DFT[x_1(n)]");

xlabel("ω/π");ylabel("幅度");

subplot(1,2,2);

stem(X1k16,""); %绘制16点 DFT的幅频特性图

title("(1b)16点DFT[x_1(n)]");

xlabel("ω/π");

ylabel("幅度");

figure(2);

subplot(2,2,1);

stem(X2k8,""); %绘制8点 DFT 的幅频特性图

title("(2a) 8点DFT[x_2(n)]");

xlabel("ω/π");

ylabel("幅度");

subplot(2,2,2);

stem(X2k16,""); %绘制16点 DFT的幅频特性图

title("(2b)16点DFT[x_2(n)]");

xlabel("ω/π");

ylabel("幅度");

subplot(2,2,3);

stem(X3k8,""); %绘制8点 DFT 的幅频特性图

title("(3a) 8点DFT[x_3(n)]");

xlabel("ω/π");

ylabel("幅度");

subplot(2,2,4);

stem(X3k16,""); %绘制16点 DFT的幅频特性图

title("(3b)16点DFT[x_3(n)]");

xlabel("ω/π");

ylabel("幅度");

运行结果图：

（2）对以下周期序列进行谱分析

πx 4(n ) =cos n ； 4

N=8;n=0:N-1; %FFT的变换区间N=8

x4n=cos(pin/4);

x5n=cos(pin/4)+cos(pin/8);

X4k8=fft(x4n); %计算x4n 的 8点DFT

X5k8=fft(x5n); %计算x5n 的 8点DFT

N=16;n=0:N-1; %FFT 的变换区间N=16

x4n=cos(pin/4);

x5n=cos(pin/4)+cos(pin/8);

X4k16=fft(x4n); %计算x4n 的 16点DFT

X5k16=fft(x5n); %计算x5n 的 16点DFT

figure(3)

subplot(2,2,1);

stem(X4k8,""); %绘制8点 DFT 的幅频特性图

title("(4a) 8点DFT[x_4(n)]");

xlabel("ω/π");

ylabel("幅度");

subplot(2,2,3);

stem(X4k16,""); %绘制16点 DFT的幅频特性图

title("(4b)16点DFT[x_4(n)]");

xlabel("ω/π");

ylabel("幅度");

subplot(2,2,2);

stem(X5k8,""); %绘制8点 DFT 的幅频特性图

title("(5a) 8点DFT[x_5(n)]");

xlabel("ω/π");

ylabel("幅度");

tisubplot(2,2,4);

stem(X5k16,""); %绘制16点 DFT的幅频特性图

title("(5b)16点DFT[x_5(n)]");

xlabel("ω/π");

ylabel("幅度");

运行结果图如下：

(3) 模拟周期信号谱分析

实验程序：

Fs=64;T=1/Fs;

N=16;n=0:N-1; %FFT 的变换区间N=16

x6nT=cos(8pinT)+cos(16pinT)+cos(20pinT); %对x6(t)16点采样 X6k16=fft(x6nT); %计算x6nT 的16点 DFT

X6k16=fftshift(X6k16); %将零频率移到频谱中心

Tp=NT;F=1/Tp; %频率分辨率F

k=-N/2:N/2-1;fk=kF; %产生16点DFT 对应的采样点频率（以零频率为中心） subplot(3,1,1);

stem(fk,abs(X6k16),"");

box on %绘制16点 DFT的幅频特性图

title("(6a) 16点|DFT[x_6(nT)]|");

xlabel("f(Hz)");

ylabel("幅度");

axis([-NF/2-1,NF/2-1,0,12max(abs(X6k16))]);

N=32;n=0:N-1; %FFT 的变换区间N=32

x6nT=cos(8pinT)+cos(16pinT)+cos(20pinT); %对x6(t)32点采样

X6k32=fft(x6nT); %计算x6nT 的32点 DFT

X6k32=fftshift(X6k32); %将零频率移到频谱中心 Tp=NT;F=1/Tp; %频率分辨率F

k=-N/2:N/2-1;fk=kF; %产生32点DFT 对应的采样点频率（以零频率为中心） subplot(3,1,2);

stem(fk,abs(X6k32),"");

box on %绘制32点 DFT的幅频特性图

title("(6b) 32点|DFT[x_6(nT)]|");

xlabel("f(Hz)");

ylabel("幅度");

axis([-NF/2-1,NF/2-1,0,12max(abs(X6k32))])

N=64;n=0:N-1; %FFT 的变换区间N=64

x6nT=cos(8pinT)+cos(16pinT)+cos(20pinT); %对x6(t)64点采样 X6k64=fft(x6nT); %计算x6nT 的64点 DFT

X6k64=fftshift(X6k64); %将零频率移到频谱中心

Tp=NT;F=1/Tp; %频率分辨率F

k=-N/2:N/2-1;fk=kF; %产生64点DFT 对应的采样点频率（以零频率为中心） subplot(3,1,3);

stem(fk,abs(X6k64),"");

box on%绘制64点 DFT的幅频特性图

title("(6a) 64点|DFT[x_6(nT)]|");

xlabel("f(Hz)");

ylabel("幅度");

axis([-NF/2-1,NF/2-1,0,12max(abs(X6k64))])

运行结果图：

实验四 IIR 数字滤波器设计

1． 实验目的

（1）熟悉用双线性变换法设计 IIR 数字滤波器的原理与方法；

（2）学会调用 MATLAB 信号处理工具箱中滤波器设计函数（或滤波器设计分析工具

fdatool ）设计各种 IIR 数字滤波器，学会根据滤波需求确定滤波器指标参数。

（3）掌握 IIR 数字滤波器的 MATLAB 实现方法。

（3）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。

2． 实验内容

（ 1）调用 buttord 和 butter 函数设计模拟低通巴特沃斯低通滤波器。 设通带截止频率 f c = 5kHz ，允许的最大衰减α p = 2dB ，阻带边缘频率 f s =12kHz ，允许的最小衰

减αs = 30dB ，试设计模拟低通巴特沃斯低通滤波器。（教材 P159，例 621） 实验程序

clc

% % 1)、给出技术指标

fp=5000; rp=1; fs=12000;rs=30;

wp=2pifp;

ws=2pifs;

% % 2)、确定巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc

% %直接计算求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc

% k_sp=sqrt((10^(01rs)-1)/(10^(01rp)-1))%求 k_sp

% labmda_sp=ws/wp %求 labmda_sp

% N=ceil(log10(k_sp)/log10(labmda_sp)) %求阶数 N

% wc=wp((10^(01rp)-1)^(-1/(2N))) %求 3dB 截止频率 wc

% % %或者调用 buttord 函数求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc

[N,wc]=buttord(wp,ws, rp, rs,"s")

%调用 buttord 函数直接求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc % %%3)、调用 buttord 函数直接求解实际的系统函数的系数

[b,a]=butter(N,wc,"s")

运行结果：

N =

5

wc =

37792e+004

b =

10e+022

0 0 0 0 0 77094

a =

10e+022

00000 00000 00000 00000 00007 77094

（ 2）用脉冲响应不变法和双线性法设计数字低通滤波器

设计低通数字滤波器，要求通带频率低于 02π rad ，允许的最大衰减 αp = 1dB ，在频

率 03π 到π 之间的阻带频率衰减大于15dB ，采样间隔 T=1s, 试用脉冲响应不变法和双线性法设计数字滤波器。（教材 P187，例 642）

实验程序：

clc

% % 1、根据数字技术指标求出模拟技术指标

T=1;

Fs=1/T;

wp=02pi;

ws=03pi;

rp=1;rs=15;

Omegap=(2/T)tan(wp/2);

%根据用双线性变换法模拟频率与数字频率之间的关系 Omega=(2/T)tan(w/2) Omegas=(2/T)tan(ws/2);

Omegap1=(wp/T);

%根据用脉冲响应不变换法模拟频率与数字频率之间的关系 Omega=(w/T) Omegas1=(ws/T);

% % 2、确定巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc

% % %调用 buttord 函数求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc

[N,wc]=buttord(Omegap,Omegas, rp, rs,"s")

%调用 buttord 函数直接求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc

[N1,wc1]=buttord(Omegap1,Omegas1, rp, rs,"s")

%调用 buttord 函数直接求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc % % 3、求实际的系统函数

[b,a]=butter(N,wc,"s")

[b1,a1]=butter(N1,wc1,"s")

% % 4、用脉冲响应不变换法和双线性变换法将模拟系统函数转化为数字系统函数

[bz,az]=bilinear(b,a,Fs)

[bz1,az1]=impinvar(b,a,Fs)

% % 5、作出该数字滤波器的幅度和相位图形

[H,f]=freqz(bz,az,256,Fs)

[H1,f]=freqz(bz1,az1,256,Fs)

% f = w/(2pi)Fs;

plot(f,20log10(abs(H)),"-",f,20log10(abs(H1)),"-");

legend("脉冲响应不变法"," 双线性法")

axis([0,03,-80,10]);

grid on;

xlabel("频率/Hz ")

ylabel("幅值/dB")

实验五 FIR 滤波器设计

一、实验目的

1、掌握 MATLAB 软件的界面使用和各种功能。

2、掌握利用 MATLAB 进行 FIR 低通滤波器的设计

3、对含噪信号进行滤波

4、了解各种窗函数对滤波器特性的影响

二、实验内容

1、用窗函数法设计数字 FIR 低通滤波器，满足下列指标： 通带截止频率：wp=02pi，通带最大衰减：Ap=025db 阻带截止频率: ws=03pi，阻带最大衰减：As=50db

并通过一信号序列的验证。

实验程序如下：

wp=02pi;wr=03pi;

wid=wr-wp;

N=ceil(66pi/wid)+1

n=[0:1:N-1];

wc=(wr+wp)/2

alpha=(N-1)/2;

m=n-alpha+eps;

hd=sin(wcm)/(pim);

ham=(hamming(N))";

hn=hdham;

[H,w]=freqz(hn,[1],1000,"whole"); Hn=(H(1:501))";wn=(w(1:501))"; mag=abs(Hn);

db=20log10((mag+eps)/max(mag)); pha=angle(Hn);

del=2pi/1000;

Rp=-(min(db(1:1:wp/del+1)))

As=-round(max(db((wr/del+1):1:501))) figure(1)

subplot(2,2,1);

plot(n,hd);

title("理想脉冲响应");

subplot(2,2,2);

plot(n,ham);

title("哈明窗");

subplot(2,2,3);

plot(n,hn);

title("实际脉冲响应");

subplot(2,2,4);

plot(wn/pi,db);

title("幅度响应(单位：DB)");

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

t=0:1/600:1;

x1=sin(2pi15t);

x2=05sin(2pi90t);

x3=02sin(2pi300t);

sig=x1+x2+x3;

newsig=fftfilt(hn,sig);

figure(2)

subplot(2,1,1);

plot(sig);

title("mixed signal");

subplot(2,1,2);

plot(newsig);

title("filted signal");

运行结果图如下：

2、用矩形窗、汉宁窗和布莱克曼窗设计 FIR 低通滤波器：

设 N=11，wc=02pi rad

实验程序如下：

N=11;Wc=02pi;

%n=0:N-1;a=(N-1)/2;

%Hdn=sin(Wc(n-a))/pi/(n-a);

%if rem(N,2)~=0

% Hdn(a+1)=Wc/pi;

%end

%wn1=boxcar(N);

%hn1=Hdnwn1"

hn1=fir1(N-1,Wc/pi,boxcar(N))

%wn2=hamming(N);

%hn2=Hdnwn2"

hn2=fir1(N-1,Wc/pi,hamming(N))

%wn3=blackman(N);

%hn3=Hdnwn3"

hn3=fir1(N-1,Wc/pi,blackman(N))

%wn4=hanning(N);

%hn4=Hdnwn4"

hn4=fir1(N-1,Wc/pi,hanning(N))

hw1=fft(hn1,1024);

hw2=fft(hn2,1024);

hw3=fft(hn3,1024);

hw4=fft(hn4,1024);

w=2[0:1023]/1024;

figure

plot(w,20log10(abs(hw1)))

hold on

plot(w,20log10(abs(hw2)) ,"-r")

plot(w,20log10(abs(hw3)) ,":g")

plot(w,20log10(abs(hw4)) ,"--m","Linewidth",2)

legend("矩形窗"," 哈明窗"," 布莱克曼窗"," 汉宁窗",1)

xlabel("Frequency in rad/sample")

ylabel("Magnitude in dB")

title("用窗函数法设计 FIR 滤波器的低通幅度特性")

运行结果图如下：

利用FFT 和IFFT 计算线性卷积

一、实验目的

学习用FFT 和IFFT 计算线性卷积的方法，并编写MATLAB 程序实现线性卷积。

二、实验内容

利用MATLAB 软件调用FFT 和IFFT 函数，编制FFT 计算线性卷积程序，并比较离散卷积计算速度。

（1）、计算线性卷积的程序

x=[1,1,1,1];

[m1,n1]=size(x);

h=[1,1,1,1,1,1];

[m2,n2]=size(h);

m=n1+n2-1;

XK=fft(x,m);

HK=fft(h,m);

YK=XKHK;

y=ifft(YK)

yc=conv(x,h)

subplot(2,1,1);

stem(real(y));title("利用FFT 进行卷积计算");

subplot(2,1,2);

stem(yc);title("利用CONV 进行卷积计算");

实验结果：

y =

10000 20000 30000 40000 40000 40000 30000 20000 10000

yc =

1 2 3 4 4 4 3 2 1

运行结果图如下：

2、比较利用函数conv 与FFT 进行线性卷积的速度

计算长度为1000的序列，x1=05n,x2=2n

（1）、实验程序

L=1000;

N=L2-1;

n=1:L;

x1=05n;

x2=2n;

t0=clock;

yc=conv(x1,x2);

tc=etime(clock,t0)

t0=clock;

yf=ifft(fft(x1,N)fft(x2,N));

tf=etime(clock,t0)

n1=0:length(yf)-1;

subplot(2,1,1);plot(n1,yc,"r");title("利用FFT 进行卷积计算");

subplot(2,1,2);plot(n1,real(yf),"b");title("利用CONV 进行卷积计算"); 实验结果：

tc =

00040

tf =

00060

tc =

00080

tf =

00070

运行结果图如下：