数字信号处理实验

实验一：信号、系统及系统响应

一、实验目的：

(1) 熟悉连续信号经理想采样前后的频谱变化关系， 加深对时域采样定理的理解。

(2) 熟悉时域离散系统的时域特性。

(3) 利用卷积方法观察分析系统的时域特性。

(4) 掌握序列傅里叶变换的计算机实现方法， 利用序列的傅里叶变换对连续信号、 离散信号及系统响应进行频域分析。 

二、实验原理：

(1) 时域采样。

(2) LTI系统的输入输出关系。

三、 实验内容 

(1) 认真复习采样理论、 离散信号与系统、 线性卷积、 序列的傅里叶变换及性质等有关内容， 阅读本实验原理与方法。 

(2) 编制实验用主程序及相应子程序。 

① 信号产生子程序， 用于产生实验中要用到的下列信号序列： 

a. xa(t)=A*e^-at *sin(Ω0t)u(t) 

A=444.128;a=50*sqrt(2)*pi; 

b. 单位脉冲序列：xb(n)=δ(n) 

c. 矩形序列： xc(n)=RN(n), N=10 

② 系统单位脉冲响应序列产生子程序。 本实验要用到两种FIR系统。 

a. ha(n)=R10(n); 

b. hb(n)=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3) 

③ 有限长序列线性卷积子程序 

用于完成两个给定长度的序列的卷积。 可以直接调用MATLAB语言中的卷积函数conv。 conv用于两个有限长度序列的卷积， 它假定两个序列都从n=0 开始。 调用格式如下：y=conv (x, h)

四、实验步骤

调通并运行实验程序， 完成下述实验内容：

 ① 分析采样序列的特性。

 a. 取采样频率fs=1 kHz, 即T=1 ms。  

b. 改变采样频率， fs=300 Hz， 观察|X(ejω)|的变化， 并做记录(打印曲线)； 进一步降低采样频率， fs=200 Hz， 观察频谱混叠是否明显存在， 说明原因， 并记录(打印)这时的|X(ejω)|曲线。

② 时域离散信号、 系统和系统响应分析。

a. 观察信号xb(n)和系统hb(n)的时域和频域特性； 利用线性卷积求信号xb(n)通过系统hb(n)的响应y(n)， 比较所求响应y(n)和hb(n)的时域及频域特性， 注意它们之间有无差别， 绘图说明， 并用所学理论解释所得结果。

b. 观察系统ha(n)对信号xc(n)的响应特性。 

③ 卷积定理的验证。

五、思考题

  (1) 在分析理想采样序列特性的实验中， 采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同? 它们所对应的模拟频率是否相同? 为什么?

  (2) 在卷积定理验证的实验中， 如果选用不同的频域采样点数M值， 例如， 选M=10和M=20， 分别做序列的傅里叶变换， 求得的结果有无差异？为什么？

六、实验程序及结果

本实验使用自定义函数的方法产生信号：

门函数：

function[y,n]=gate(np,ns,nf)

if ns>np|ns>nf|np>nf

    error('ÊäÈëλÖòÎÊý²»Âú×ãns<=np<=nf');

else

    n=[ns:nf];

y=(n>=0&nend

冲击函数：

function[y,n]=impseq(np,ns,nf)

if ns>np|ns>nf|np>nf

    error('ÊäÈëλÖòÎÊý²»Âú×ãns<=np<=nf');

else

    n=[ns:nf];

    y=[(n-np)==0];

end

时域函数：

function[y,ts]=sig(t0,tp,t1);

A=444.128;a=50.*sqrt(2).*pi;w0=50.*sqrt(2).*pi;

ts=t0:tp:t1;

y=A.*exp(-a.*ts).*sin(w0.*ts).*(ts>=0);

;hb=impseq(0,0,10)+2.5.*impseq(1,0,10)+2.5.*impseq(2,0,10)+impseq(3,0,10);

;[ha,N]=gate(10,0,10);

第一部分主程序：

[y1,t1]=sig(0,0.001,1);

%t1=0:0.001:1;

%y1=0.5.*sin(2.*pi.*15.*t1);

figure;

plot(t1,y1);

axis([0 0.15 -10 150]);

title('f=1000ʱÓòÐźÅ');

figure;

f=-500:500;

yy1=fftshift(fft(y1));

plot(f,abs(yy1));

title('f=1000ƵÆ×');

[y2,t2]=sig(0,1/300,1);

figure;

plot(t2,y2);

axis([0 0.15 -10 150]);

title('f=300ʱÓòÐźÅ');

figure;

yy2=fftshift(fft(y2));

plot(-150:150,abs(yy2));

title('f=300ƵÆ×');

[y3,t3]=sig(0,0.02,1);

figure;

plot(t3,y3);

title('f=200ʱÓòÐźÅ');

axis([0 0.15 -8 5]);

figure;

yy3=fftshift(fft(y3));

plot(-25:25,abs(yy3));

title('f=200ƵÆ×');

第二部分主程序：

xb=impseq(0,0,100);

figure;

stem(0:length(xb)-1,xb);

axis([-1,3,0,2]);

figure;

y1=fft(xb);

plot(0:length(y1)-1,y1);

hb=impseq(0,0,100)+2.5.*impseq(1,0,100)+2.5.*impseq(2,0,100)+impseq(3,0,100);

figure;

stem(0:length(hb)-1,hb);

axis([-1,5,0,3])

figure

y2=fftshift(fft(hb));

plot(-(length(y2)-1)/2:(length(y2)-1)/2,abs(y2));

figure;

y3=conv(double(xb),double(hb));

plot(0:length(y3)-1,y3);

y4=fftshift(fft(y3));

figure;

plot(-(length(y4)-1)/2:(length(y4)-1)/2,abs(y4));

xc=gate(10,0,100);

ha=xc;

y5=fftshift(fft(conv(double(xc),double(ha))));

figure;

plot(-(length(y5)-1)/2:(length(y5)-1)/2,abs(y5));

实验结果：

第一部分：

采样频率为200Hz时时域恢复                         200Hz时频谱

采样频率为300Hz时时域恢复                         300Hz时频谱

采样频率为1000Hz是时域恢复                        1000Hz时频谱

第二部分：

Xb频域分析                                 xb时域分析

hb频域分析                                 hb时域分析

信号通过xb频域                     信号通过xc频域

卷积定理的验证：

%%%%%卷积定理的验证

y1=[0,1,2,3,4,5,6,7,8];

y2=[8,7,6,5,4,3,2,1,0];

y3=conv(y1,y2);

f1=fft(y1,20);

f2=fft(y2,20);

f3=fft(y3,20);

f=f1.*f2;

subplot(2,1,1);

stem(0:length(f)-1,f,'.');

subplot(2,1,2);

stem(0:length(f3)-1,f3,'.');

时域卷积后求频谱和频域相乘

可见，f= conv(y1,y2)的频谱和y1,y2的频谱相乘后结果相同。即满足卷积定理

七、回答问题。

(1) 在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？它们所对应的模拟频率是否相同？为什么？ 

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

(2) 在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得的结果有无差异？ 

答：有差异，所到的结果点数不同。

八、总结。

本实验主要是后续实验的基础。涉及内容也比较浅显。

单位冲击序列与hb卷积后得到的频谱与hb原频谱相同。原因很简单，是因为单位冲击序列卷积任何函数仍然是原函数。

实验二： 用FFT作谱分析

一、实验目的

(1) 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法， 所以FFT的运算结果必然满足DFT的基本性质)。

(2) 熟悉FFT算法原理和FFT子程序的应用。

(3) 学习用FFT对连续信号和时域离散信号进行谱分析的方法， 了解可能出现的分析误差及其原因， 以便在实际中正确应用FFT。

二、实验步骤

(1) 复习DFT的定义、 性质和用DFT作谱分析的有关内容。

(2) 复习FFT算法原理与编程思想， 并对照DIT-FFT运算流图和程序框图， 读懂本实验提供的FFT子程序。

  (3) 编制信号产生子程序， 产生以下典型信号供谱分析用：

(4) 编写主程序。

下图给出了主程序框图， 供参考。 本实验提供FFT子程序和通用绘图子程序。

(5) 按实验内容要求， 上机实验， 并写出实验报告。 

三、实验内容

(1) 对 2 中所给出的信号逐个进行谱分析。          

(2) 令x(n)=x4(n)+x5(n)， 用FFT计算 8 点和 16 点离散傅里叶变换，

   X(k)=DFT［x(n)］

(3) 令x(n)=x4(n)+jx5(n)， 重复(2)。

四、思考题

  (1) 在N=8时， x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢?

  (2) 如果周期信号的周期预先不知道， 如何用FFT进行谱分析?

五、实验程序及结果

x1 = [1,1,1,1,0,0,0,0,0,0];

x2 = [1,2,3,4,4,3,2,1,0,0];

x3 = [4,3,2,1,1,2,3,4,0,0];

n=0:1000;

x4=cos(n);

x5=sin(n);

t=0:0.001:1;

x6=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t);

% figure;

% title('Öð¸öÆ×·ÖÎö')

subplot(3,2,1);

stem(0:length(x1)-1,x1);

subplot(3,2,2);

stem(0:length(x2)-1,x2);

subplot(3,2,3);

stem(0:length(x3)-1,x3);

subplot(3,2,4);

stem(0:20,x4(1:21));

subplot(3,2,5);

stem(0:20,x5(1:21));

subplot(3,2,6);

stem(0:20,x6(1:21));

figure;

stem(0:200,x6(1:201));

y1=abs(fft(x1,8));

y2=abs(fft(x2,8));

y3=abs(fft(x3,8));

y4=abs(fft(x4,8));

y5=abs(fft(x5,8));

y6=abs(fft(x6,8));

figure;

subplot(3,2,1);

stem(0:length(y1)-1,y1);

subplot(3,2,2);

stem(0:length(y2)-1,y2);

subplot(3,2,3);

stem(0:length(y3)-1,y3);

subplot(3,2,4);

stem(0:length(y1)-1,y1);

subplot(3,2,5);

stem(0:length(y5)-1,y5);

subplot(3,2,6);

stem(0:length(y6)-1,y6);

%%%figure;

%%%stem(0:15,y6(1:16));

%%%%%%%%%%%%%%%

y1=abs(fft(x1,16));

y2=abs(fft(x2,16));

y3=abs(fft(x3,16));

y4=abs(fft(x4,16));

y5=abs(fft(x5,16));

y6=abs(fft(x6,16));

figure;

subplot(3,2,1);

stem(0:length(y1)-1,y1);

subplot(3,2,2);

stem(0:length(y2)-1,y2);

subplot(3,2,3);

stem(0:length(y3)-1,y3);

subplot(3,2,4);

stem(0:length(y6)-1,y4);

subplot(3,2,5);

stem(0:length(y6)-1,y5);

subplot(3,2,6);

stem(0:length(y6)-1,y6);

%figure;

%stem(0:15,y6(1:16));

x7=x4+j*x5;

y71=abs(fft(x7,8));

y72=abs(fft(x7,16));

figure;

subplot(1,2,1);

stem(0:length(y71)-1,y71);

subplot(1,2,2);

stem(0:length(y72)-1,y72);

程序结果如下所示：

每个信号时域分析

因为第六个信号不易在上图中展示（包含信息不到一个周期）所以单独画出图形如下：

x6的时域波形

逐个进行8点离散傅里叶变换

逐个进行16点离散傅里叶变换

x(n)8点和16点离散傅里叶变换比较

六、回答问题

(1) 在N=8时，x2(n)和x3(n)的幅频特性会相同吗? 为什么? N=16呢?

答：N=8时幅频特性一样，N=16时幅频特性不一样。Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。

总结：采样点N的不同，分辨率不同，在N=8，俩函数不同的部分没有分辨出来，因此特性相同，而N=16时已经足以分辨出不同的频率，因此频谱不同。

(2) 如果周期信号的周期预先不知道，如何用FFT进行谱分析?

答：设一个定长的m值,先取2m,看2m/m的误差是否大，如大的话再取4m,看4m/2m的误差是否大，如不大，4m(4倍的m值)则可近似原来点的谱分析。

七、总结

    本实验主要掌握fft函数的用法，包括采样点的意义。这个意义已经在七、回答问题中描述，不再赘述。

实验三：用窗函数法设计FIR数字滤波器

一、实验目的

(1) 掌握用窗函数法设计FIR数字滤波器的原理和方法。

(2) 熟悉线性相位FIR数字滤波器特性。

(3) 了解各种窗函数对滤波特性的影响。

二、实验内容及步骤

  (1) 复习用窗函数法设计FIR数字滤波器一节内容， 阅读本实验原理， 掌握设计步骤。

  (2) 编写程序。

  ① 编写能产生矩型窗、 升余弦窗、 改进升余弦窗和二阶升余弦窗的窗函数子程序。

  ② 编写主程序。 其中幅度特性要求用dB表示。

窗函数法设计滤波器主程序框图

三、思考题

  (1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减， 如何用窗函数法设计线性相位低通滤波器? 写出设计步骤。

  (2) 如果要求用窗函数法设计带通滤波器， 且给定上、 下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应hd(n)。 

四、实验程序及结果

myfreqz(b,a)函数，是自编的改进freqz()的函数。

function [db,mag,pha,grd,w] =myfreqz(b,a);

%advanced freqz

%[db,mag,pha,grd,w] =myfreqz(b,a);

[H,w]=freqz(b,a,1000,'whole');

H=(H(1:1:501))';w=(w(1:1:501))';

mag=abs(H);

db=20*log10((mag+eps)/max(mag));

pha=angle(H);

grd=grpdelay(b,a,w);

以下是主程序部分

wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp;

N0=ceil(6.6*pi/deltaw);

N=N0+mod(N0+1,2);

wdham=(hamming(N))';

wc=(ws+wp)/2;

% hd=ideallp(wc,N);

tao=(N-1)/2;

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

m=n-tao+eps;

hd=sin(wc*m)./(pi*m);

h=hd.*wdham;

[db,mag,pha,grd,w]=myfreqz(h,[1]);

dw=2*pi/1000;

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

As=-round(max(db(ws/dw+1:501)));

figure;

subplot(2,2,1);

stem(0:length(hd)-1,hd,'.');

subplot(2,2,2);

stem(0:length(wdham)-1,wdham,'.');

subplot(2,2,3);

stem(0:length(h)-1,h,'.');

subplot(2,2,4);

plot(1/length(db):1/length(db):1,db);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

wn1=boxcar(N);

figure;

subplot(2,2,1);

stem(0:N-1,wn1,'.');

title('boxcar');

wn2=bartlett(N);

subplot(2,2,2);

stem(0:N-1,wn2,'.');title('bartlett');

wn3=hanning(N);

subplot(2,2,3);

stem(0:N-1,wn3,'.');title('hanning');

% wn4=hamming(N);

wn4=blackman(N);

subplot(2,2,4);

stem(0:N-1,wn4,'.');title('blackman');

% wn=(N,beta);

得到的结果如下图：

其中，图一是理想脉冲响应，第一行图二是汉明窗；图三是实际脉冲响应，图四是幅度响应。一下是几个典型窗函数的图像。，汉明窗上例已画出，不再画。每个窗的名字已在每个子图上方标出。

典型窗函数

五、回答问题

(1) 给定通带截止频率和阻带截止频率以及阻带最小衰减，用窗函数法设计线性相位低通滤波器的设计步骤： 

答： 

技术指标 Wp=0.2*pi,Ws=0.4*pi,Ap=0.25dB,As=50dB

方法一

选择海明窗

clear all;

Wp=0.2*pi;

Ws=0.4*pi;

tr_wide=Ws-Wp;    %过渡带宽度

N=ceil(6.6*pi/tr_wide)+1;   %滤波器长度

n=0:1:N-1;

Wc=(Wp+Ws)/2;                %理想低通滤波器的截止频率

hd=ideal_lp1(Wc,N);            %理想滤波器的单位冲击响应

w_ham=(hamming(N))';      %海明窗

h=hd.*w_ham;              %实际海明窗的响应

[db,mag,pha,w]=freqz_m2(h,[1]);  %计算实际滤波器的幅度响应

delta_w=2*pi/1000;

Ap=-(min(db(1:1:Wp/delta_w+1)))   %实际通带纹波

As=-round(max(db(Ws/delta_w+1:1:501)))    %实际阻带纹波

subplot(221)

stem(n,hd)

title('理想单位脉冲响应hd(n)')

subplot(222)

stem(n,w_ham)

title('海明窗')

subplot(223)

stem(n,h)

title('实际单位脉冲响应hd(n)')

subplot(224)

plot(wi/pi,db)

title('幅度响应(dB)')

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

方法二

Window=blackman(16);

b=fir1(15,0.3*pi ,'low',Window);

freqz(b,128)

 (2)用窗函数法设计带通滤波器， 且给定上、 下边带截止频率为ω1和ω2，理想带通的单位脉冲响应hd(n)的求解过程： 

答：由hd(n)=1/2π∫w1woe-jwaejnwdw+1/2π∫w2w3e-jwaejnwdw 

（其中w0=-w0-wc,w1=-w0+wc,w2=w0-wc,w3=w0+wc） 

计算整理后可得： 

hd（n）=2/((n-a)*π)*sin[(n-a)wc]*cos[(n-a)w0] 

=2wc/π* sa[(n-a)wc]*cos[(n-a)w0]