实验一 系统响应及系统稳定性
一、实验目的
（1）掌握 求系统响应的方法。
（2）掌握时域离散系统的时域特性。
（3）分析、观察及检验系统的稳定性。
二、实验内容
（1）给定一个低通滤波器的差分方程为
y(n)=0.05x(n)+0.05x(n-1)+0.9y(n-1)
输入信号x1(n)=R8(n)
x2(n)=u(n)
(a) 分别求出系统对x1(n)=R8(n) 和x2(n)=u(n)的响应序列，并画出其波形。
(b) 求出系统的单位冲响应，画出其波形。
实验程序：
A=[1,-0.9];B=[0.05,0.05]; %%系统差分方程系数向量 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)+2.5δ(n-1)+2.5δ(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 2.5 2.5 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)=1.8237y(n-1)-0.9801y(n-2)+b0x(n)-b0x(n-2)
令b0=1/100.49，谐振器的谐振频率为 0.4rad。
(a) 用实验方法检查系统是否稳定。输入信号为u(n) 时，画出系统输出波形。
(b) 给定输入信号为
x(n)= sin(0.014n )+ sin(0.4n )
求出系统的输出响应，并画出其波形。
实验程序：
un=ones(1,256); %产生信号 u(n)
n=0:255;
xsin=sin(0.014*n)+sin(0.4*n); %产生正弦信号
A=[1,-1.8237,0.9801];
B=[1/100.49,0,-1/100.49]; %系统差分方程系数向量 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=444.128，α =50 2 π，Ω =50 2 πrad/s
Tp=64/1000; %观察时间 Tp=64 微秒
Fs=1000;
T=1/Fs;
M=Tp*Fs;
n=0:M-1;
t=n*T;
A=444.128;
alph=pi*50*2^0.5;
omega=pi*50*2^0.5;
xat=A*exp(-alph.*t).*sin(omega*t);%给定模拟信号
Xk=T*fft(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) T*FT[x_a(nT)],F_s=1000Hz');
xlabel('\omega/hz');
ylabel('(H_1(e^j^w))');
axis([0,Fs,0,1.2*max(abs(Xk))]);
Fs=300;T=1/Fs;
M=Tp*Fs;
n=0:M-1;
t=n*T;
A=444.128;
alph=pi*50*2^0.5;
omega=pi*50*2^0.5;
xat=A*exp(-alph*t).*sin(omega*t);
Xk=T*fft(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) T*FT[x_a(nT)],Fs=300Hz');
xlabel('\omega/hz');
ylabel('(H_2(e^j^w))');
axis([0,Fs,0,1.2*max(abs(Xk))]);
A=444.128;
alph=pi*50*2^0.5;
omega=pi*50*2^0.5;
xat=A*exp(-alph*t).*sin(omega*t);
Xk=T*fft(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) T*FT[x_a(nT)],Fs=200Hz');
xlabel('\omega/hz');
ylabel('(H_3(e^j^w))');
axis([0,Fs,0,1.2*max(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=2*k/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(pi*n/4);
x5n=cos(pi*n/4)+cos(pi*n/8);
X4k8=fft(x4n); %计算x4n 的 8点DFT
X5k8=fft(x5n); %计算x5n 的 8点DFT
N=16;n=0:N-1; %FFT 的变换区间N=16
x4n=cos(pi*n/4);
x5n=cos(pi*n/4)+cos(pi*n/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(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)16点采样 X6k16=fft(x6nT); %计算x6nT 的16点 DFT
X6k16=fftshift(X6k16); %将零频率移到频谱中心
Tp=N*T;F=1/Tp; %频率分辨率F
k=-N/2:N/2-1;fk=k*F; %产生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([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))]);
N=32;n=0:N-1; %FFT 的变换区间N=32
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)32点采样
X6k32=fft(x6nT); %计算x6nT 的32点 DFT
X6k32=fftshift(X6k32); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率F
k=-N/2:N/2-1;fk=k*F; %产生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([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))])
N=64;n=0:N-1; %FFT 的变换区间N=64
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)64点采样 X6k64=fft(x6nT); %计算x6nT 的64点 DFT
X6k64=fftshift(X6k64); %将零频率移到频谱中心
Tp=N*T;F=1/Tp; %频率分辨率F
k=-N/2:N/2-1;fk=k*F; %产生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([-N*F/2-1,N*F/2-1,0,1.2*max(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，例 6.2.1） 实验程序
clc
% % 1)、给出技术指标
fp=5000; rp=1; fs=12000;rs=30;
wp=2*pi*fp;
ws=2*pi*fs;
% % 2)、确定巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc
% %直接计算求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc
% k_sp=sqrt((10.^(0.1.*rs)-1)/(10.^(0.1.*rp)-1))%求 k_sp
% labmda_sp=ws/wp %求 labmda_sp
% N=ceil(log10(k_sp)/log10(labmda_sp)) %求阶数 N
% wc=wp*((10.^(0.1.*rp)-1).^(-1/(2*N))) %求 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 =
3.7792e+004
b =
1.0e+022 *
0 0 0 0 0 7.7094
a =
1.0e+022 *
0.0000 0.0000 0.0000 0.0000 0.0007 7.7094
（ 2）用脉冲响应不变法和双线性法设计数字低通滤波器
设计低通数字滤波器，要求通带频率低于 0.2π rad ，允许的最大衰减 αp = 1dB ，在频
率 0.3π 到π 之间的阻带频率衰减大于15dB ，采样间隔 T=1s, 试用脉冲响应不变法和双线性法设计数字滤波器。（教材 P187，例 6.4.2）
实验程序：
clc
% % 1、根据数字技术指标求出模拟技术指标
T=1;
Fs=1/T;
wp=0.2*pi;
ws=0.3*pi;
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/(2*pi)*Fs;
plot(f,20*log10(abs(H)),'-.',f,20*log10(abs(H1)),'-');
legend('脉冲响应不变法',' 双线性法')
axis([0,0.3,-80,10]);
grid on;
xlabel('频率/Hz ')
ylabel('幅值/dB')
实验五 FIR 滤波器设计
一、实验目的
1、掌握 MATLAB 软件的界面使用和各种功能。
2、掌握利用 MATLAB 进行 FIR 低通滤波器的设计
3、对含噪信号进行滤波
4、了解各种窗函数对滤波器特性的影响
二、实验内容
1、用窗函数法设计数字 FIR 低通滤波器，满足下列指标： 通带截止频率：wp=0.2*pi，通带最大衰减：Ap=0.25db 阻带截止频率: ws=0.3*pi，阻带最大衰减：As=50db
并通过一信号序列的验证。
实验程序如下：
wp=0.2*pi;wr=0.3*pi;
wid=wr-wp;
N=ceil(6.6*pi/wid)+1
n=[0:1:N-1];
wc=(wr+wp)/2
alpha=(N-1)/2;
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
ham=(hamming(N))';
hn=hd.*ham;
[H,w]=freqz(hn,[1],1000,'whole'); Hn=(H(1:501))';wn=(w(1:501))'; mag=abs(Hn);
db=20*log10((mag+eps)/max(mag)); pha=angle(Hn);
del=2*pi/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(2*pi*15*t);
x2=0.5*sin(2*pi*90*t);
x3=0.2*sin(2*pi*300*t);
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=0.2*pi rad
实验程序如下：
N=11;Wc=0.2*pi;
%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=Hdn.*wn1'
hn1=fir1(N-1,Wc/pi,boxcar(N))
%wn2=hamming(N);
%hn2=Hdn.*wn2'
hn2=fir1(N-1,Wc/pi,hamming(N))
%wn3=blackman(N);
%hn3=Hdn.*wn3'
hn3=fir1(N-1,Wc/pi,blackman(N))
%wn4=hanning(N);
%hn4=Hdn.*wn4'
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,20*log10(abs(hw1)))
hold on
plot(w,20*log10(abs(hw2)) ,'-.r')
plot(w,20*log10(abs(hw3)) ,':.g')
plot(w,20*log10(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=XK.*HK;
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 =
1.0000 2.0000 3.0000 4.0000 4.0000 4.0000 3.0000 2.0000 1.0000
yc =
1 2 3 4 4 4 3 2 1
运行结果图如下：
2、比较利用函数conv 与FFT 进行线性卷积的速度
计算长度为1000的序列，x1=0.5*n,x2=2*n
（1）、实验程序
L=1000;
N=L*2-1;
n=1:L;
x1=0.5*n;
x2=2*n;
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 =
0.0040
tf =
0.0060
tc =
0.0080
tf =
0.0070
运行结果图如下：
实验一 系统响应及系统稳定性
一、实验目的
（1）掌握 求系统响应的方法。
（2）掌握时域离散系统的时域特性。
（3）分析、观察及检验系统的稳定性。
二、实验内容
（1）给定一个低通滤波器的差分方程为
y(n)=0.05x(n)+0.05x(n-1)+0.9y(n-1)
输入信号x1(n)=R8(n)
x2(n)=u(n)
(a) 分别求出系统对x1(n)=R8(n) 和x2(n)=u(n)的响应序列，并画出其波形。
(b) 求出系统的单位冲响应，画出其波形。
实验程序：
A=[1,-0.9];B=[0.05,0.05]; %%系统差分方程系数向量 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)+2.5δ(n-1)+2.5δ(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 2.5 2.5 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)=1.8237y(n-1)-0.9801y(n-2)+b0x(n)-b0x(n-2)
令b0=1/100.49，谐振器的谐振频率为 0.4rad。
(a) 用实验方法检查系统是否稳定。输入信号为u(n) 时，画出系统输出波形。
(b) 给定输入信号为
x(n)= sin(0.014n )+ sin(0.4n )
求出系统的输出响应，并画出其波形。
实验程序：
un=ones(1,256); %产生信号 u(n)
n=0:255;
xsin=sin(0.014*n)+sin(0.4*n); %产生正弦信号
A=[1,-1.8237,0.9801];
B=[1/100.49,0,-1/100.49]; %系统差分方程系数向量 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=444.128，α =50 2 π，Ω =50 2 πrad/s
Tp=64/1000; %观察时间 Tp=64 微秒
Fs=1000;
T=1/Fs;
M=Tp*Fs;
n=0:M-1;
t=n*T;
A=444.128;
alph=pi*50*2^0.5;
omega=pi*50*2^0.5;
xat=A*exp(-alph.*t).*sin(omega*t);%给定模拟信号
Xk=T*fft(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) T*FT[x_a(nT)],F_s=1000Hz');
xlabel('\omega/hz');
ylabel('(H_1(e^j^w))');
axis([0,Fs,0,1.2*max(abs(Xk))]);
Fs=300;T=1/Fs;
M=Tp*Fs;
n=0:M-1;
t=n*T;
A=444.128;
alph=pi*50*2^0.5;
omega=pi*50*2^0.5;
xat=A*exp(-alph*t).*sin(omega*t);
Xk=T*fft(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) T*FT[x_a(nT)],Fs=300Hz');
xlabel('\omega/hz');
ylabel('(H_2(e^j^w))');
axis([0,Fs,0,1.2*max(abs(Xk))]);
A=444.128;
alph=pi*50*2^0.5;
omega=pi*50*2^0.5;
xat=A*exp(-alph*t).*sin(omega*t);
Xk=T*fft(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) T*FT[x_a(nT)],Fs=200Hz');
xlabel('\omega/hz');
ylabel('(H_3(e^j^w))');
axis([0,Fs,0,1.2*max(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=2*k/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(pi*n/4);
x5n=cos(pi*n/4)+cos(pi*n/8);
X4k8=fft(x4n); %计算x4n 的 8点DFT
X5k8=fft(x5n); %计算x5n 的 8点DFT
N=16;n=0:N-1; %FFT 的变换区间N=16
x4n=cos(pi*n/4);
x5n=cos(pi*n/4)+cos(pi*n/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(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)16点采样 X6k16=fft(x6nT); %计算x6nT 的16点 DFT
X6k16=fftshift(X6k16); %将零频率移到频谱中心
Tp=N*T;F=1/Tp; %频率分辨率F
k=-N/2:N/2-1;fk=k*F; %产生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([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))]);
N=32;n=0:N-1; %FFT 的变换区间N=32
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)32点采样
X6k32=fft(x6nT); %计算x6nT 的32点 DFT
X6k32=fftshift(X6k32); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率F
k=-N/2:N/2-1;fk=k*F; %产生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([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))])
N=64;n=0:N-1; %FFT 的变换区间N=64
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)64点采样 X6k64=fft(x6nT); %计算x6nT 的64点 DFT
X6k64=fftshift(X6k64); %将零频率移到频谱中心
Tp=N*T;F=1/Tp; %频率分辨率F
k=-N/2:N/2-1;fk=k*F; %产生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([-N*F/2-1,N*F/2-1,0,1.2*max(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，例 6.2.1） 实验程序
clc
% % 1)、给出技术指标
fp=5000; rp=1; fs=12000;rs=30;
wp=2*pi*fp;
ws=2*pi*fs;
% % 2)、确定巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc
% %直接计算求巴特沃斯低通滤波器的阶数 N 以及 3dB 截止频率 wc
% k_sp=sqrt((10.^(0.1.*rs)-1)/(10.^(0.1.*rp)-1))%求 k_sp
% labmda_sp=ws/wp %求 labmda_sp
% N=ceil(log10(k_sp)/log10(labmda_sp)) %求阶数 N
% wc=wp*((10.^(0.1.*rp)-1).^(-1/(2*N))) %求 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 =
3.7792e+004
b =
1.0e+022 *
0 0 0 0 0 7.7094
a =
1.0e+022 *
0.0000 0.0000 0.0000 0.0000 0.0007 7.7094
（ 2）用脉冲响应不变法和双线性法设计数字低通滤波器
设计低通数字滤波器，要求通带频率低于 0.2π rad ，允许的最大衰减 αp = 1dB ，在频
率 0.3π 到π 之间的阻带频率衰减大于15dB ，采样间隔 T=1s, 试用脉冲响应不变法和双线性法设计数字滤波器。（教材 P187，例 6.4.2）
实验程序：
clc
% % 1、根据数字技术指标求出模拟技术指标
T=1;
Fs=1/T;
wp=0.2*pi;
ws=0.3*pi;
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/(2*pi)*Fs;
plot(f,20*log10(abs(H)),'-.',f,20*log10(abs(H1)),'-');
legend('脉冲响应不变法',' 双线性法')
axis([0,0.3,-80,10]);
grid on;
xlabel('频率/Hz ')
ylabel('幅值/dB')
实验五 FIR 滤波器设计
一、实验目的
1、掌握 MATLAB 软件的界面使用和各种功能。
2、掌握利用 MATLAB 进行 FIR 低通滤波器的设计
3、对含噪信号进行滤波
4、了解各种窗函数对滤波器特性的影响
二、实验内容
1、用窗函数法设计数字 FIR 低通滤波器，满足下列指标： 通带截止频率：wp=0.2*pi，通带最大衰减：Ap=0.25db 阻带截止频率: ws=0.3*pi，阻带最大衰减：As=50db
并通过一信号序列的验证。
实验程序如下：
wp=0.2*pi;wr=0.3*pi;
wid=wr-wp;
N=ceil(6.6*pi/wid)+1
n=[0:1:N-1];
wc=(wr+wp)/2
alpha=(N-1)/2;
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);
ham=(hamming(N))';
hn=hd.*ham;
[H,w]=freqz(hn,[1],1000,'whole'); Hn=(H(1:501))';wn=(w(1:501))'; mag=abs(Hn);
db=20*log10((mag+eps)/max(mag)); pha=angle(Hn);
del=2*pi/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(2*pi*15*t);
x2=0.5*sin(2*pi*90*t);
x3=0.2*sin(2*pi*300*t);
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=0.2*pi rad
实验程序如下：
N=11;Wc=0.2*pi;
%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=Hdn.*wn1'
hn1=fir1(N-1,Wc/pi,boxcar(N))
%wn2=hamming(N);
%hn2=Hdn.*wn2'
hn2=fir1(N-1,Wc/pi,hamming(N))
%wn3=blackman(N);
%hn3=Hdn.*wn3'
hn3=fir1(N-1,Wc/pi,blackman(N))
%wn4=hanning(N);
%hn4=Hdn.*wn4'
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,20*log10(abs(hw1)))
hold on
plot(w,20*log10(abs(hw2)) ,'-.r')
plot(w,20*log10(abs(hw3)) ,':.g')
plot(w,20*log10(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=XK.*HK;
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 =
1.0000 2.0000 3.0000 4.0000 4.0000 4.0000 3.0000 2.0000 1.0000
yc =
1 2 3 4 4 4 3 2 1
运行结果图如下：
2、比较利用函数conv 与FFT 进行线性卷积的速度
计算长度为1000的序列，x1=0.5*n,x2=2*n
（1）、实验程序
L=1000;
N=L*2-1;
n=1:L;
x1=0.5*n;
x2=2*n;
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 =
0.0040
tf =
0.0060
tc =
0.0080
tf =
0.0070
运行结果图如下：