离散时间系统及离散卷积

实验一、离散时间系统及离散卷积

1、单位脉冲响应

源程序：

function pr1() %定义函数pr1

a=[1,-1,0.9]; %定义差分方程y(n)-y(n-1)+0.9y(n-2)=x(n) b=1;

x=impseq(0,-20,120); %调用impseq函数（matlab软件的函数库）

n=[-20:120]; %定义n的范围，从-20 到120

h=filter(b,a,x); %调用函数给纵坐标赋值

figure(1) %绘图figure 1 (冲激响应)

stem(n,h); %在图中绘出冲激

title('单位冲激响应（耿海锋）'); %定义标题为:'冲激响应（耿海锋）'

xlabel('n'); %绘图横座标为n

ylabel('h(n)'); %绘图纵座标为h(n)

figure(2) %绘图figure 2

[z,p,g]=tf2zp(b,a); %绘出零极点图

zplane(z,p)

function [x,n]=impseq(n0,n1,n2) %声明impseq函数

n=[n1:n2];

x=[(n-n0)==0];

结果：

Figure 1:

Figure 2:

2、离散系统的幅频、相频的分析

源程序：

function pr2()

b=[0.0181,0.0543,0.0543,0.0181];

a=[1.000,-1.76,1.1829,-0.2781];

m=0:length(b)-1; % m的范围，从0 到3

l=0:length(a)-1; % l的范围，从0 到3

K=5000;

k=1:K;

w=pi*k/K; %角频率w

H=(b*exp(-j*m'*w))./(a*exp(-j*l'*w));%对系统函数的定义

figure(1)

magH=abs(H); %magH为幅度

angH=angle(H); %angH为相位

plot(w/pi,magH-耿海锋); %绘制w(pi)-magH-耿海锋的图形

figure(2)

axis([0,1,0,1]); %限制横纵座标从0到1

xlabel('w(pi)'); %x座标为 w(pi)

ylabel('|H|'); %y座标为 angle(H)-耿海锋

title('幅度，相位响应（耿海锋）'); %图的标题为:'幅度，相位响应（耿海锋）'

plot(w/pi,angH); %绘制w(pi)-angH的图形

grid; %为座标添加名称

xlabel('w(pi)'); %x座标为 w(pi)

ylabel('angle(H)'); %y座标为 angle(H)

结果：

Figure1

Figure2

3、卷积计算

源程序：

function pr3()

n=-5:50; %声明n的范围，从-5到50

u1=stepseq(0,-5,50); %调用stepseq函数声用明u1=u(n)

u2=stepseq(10,-5,50); %调用stepseq函数声用明u2=u(n-10)

%输入x(n)和冲激响应h(n)

x=u1-u2; %x(n)=u(n)-u(n-10)

h=((0.9).^n).*u1; %h(n)=0.9^n*u(n)

figure(1)

subplot(3,1,1); %绘制第一个子图

stem(n,x); %绘制图中的冲激

axis([-5,50,0,2]); %限定横纵座标的范围

title('输入序列-52101702-耿海锋'); %规定标题为:'输入序列-52101702-耿海锋'

xlabel('n'); %横轴为n

ylabel('x(n)'); %纵轴为x(n)

subplot(3,1,2); %绘制第二个子图

stem(n,h); %绘制图中的冲激

axis([-5,50,0,2]); %限定横纵座标的范围

title('冲激响应序列-52101702-耿海锋'); %规定标题为:'冲激响应序列-52101702-耿海锋'

xlabel('n'); %横轴为n

ylabel('h(n)'); %纵轴为h(n)

%输出响应

[y,ny]=conv_m(x,n,h,n); %调用conv_m函数

subplot(3,1,3); %绘制第三个子图

stem(ny,y);

axis([-5,50,0,8]);

title('输出响应-52101702-耿海锋'); %规定标题为:'输出响应-52101702-耿海锋'

xlabel('n');

ylabel('y(n)'); %纵轴为y(n)

%stepseq.m子程序

%实现当n>=n0时x(n)的值为1

function [x,n]=stepseq(n0,n1,n2)

n=n1:n2;

x=[(n-n0)>=0];

%con_m的子程序

%实现卷积的计算

function [y,ny]=conv_m(x,nx,h,nh) nyb=nx(1)+nh(1);

nye=nx(length(x))+nh(length(h)); ny=[nyb:nye];

y=conv(x,h);

结果：

实验二、离散傅立叶变换与快速傅立叶变换

1、离散傅立叶变换（DFT）

源程序：

function pr4()

F=50;

N=64;

T=0.000625;

n=1:N;

x=cos(2*pi*F*n*T); %x(n)=cos(pi*n/16)

subplot(2,1,1); %绘制第一个子图x(n)

stem(n,x); %绘制冲激

title('x(n)'); %标题为x(n)

xlabel('n'); %横座标为n

X=dft(x,N); %调用dft函数计算x(n)的傅里叶变换magX=abs(X); %取变换的幅值

subplot(2,1,2); %绘制第二个子图DFT|X|

stem(n,X);

title('DFT|X|');

xlabel('f(pi)'); %横座标为f(pi)

%dft的子程序

%实现离散傅里叶变换

function [Xk]=dft(xn,N)

n=0:N-1;

k=0:N-1;

WN=exp(-j*2*pi/N);

nk=n'*k;

WNnk=WN.^nk;

Xk=xn*WNnk;

结果：

F=50,N=64,T=0.000625时的波形

F=50,N=32,T=0.000625时的波形:

2、快速傅立叶变换（FFT）

源程序：

%function pr5()

F=50;N=64;T=0.000625;

n=1:N;

x=cos(2*pi*F*n*T); %x(n)=cos(pi*n/16) subplot(2,1,1);plot(n,x);

title('x(n)');xlabel('n'); %在第一个子窗中绘图x(n)

X=fft(x);magX=abs(X);

subplot(2,1,2);plot(n,X);

title('DTFT|X|');xlabel('f(pi)'); %在第二个子图中绘图x(n)的快速

傅%里叶变换

结果：

3、卷积的快速算法

源程序：

function pr6()

n=0:15;

x=1.^n;

h=(4/5).^n;

x(16:32)=0;

h(16:32)=0;

%到此 x(n)=1, n=0~15; x(n)=0,n=16~32 % h(n)=(4/5)^n, n=0~15; h(n)=0,n=16~32

subplot(3,1,1);

stem(x);

title('x(n)');

axis([1,32,0,1.5]); %在第一个子窗绘图x(n)横轴从1到32,纵轴从0到1.5（边框范围）

subplot(3,1,2);

stem(h);

title('h(n)');

axis([1,32,0,1.5]); %在第二个子窗绘图h(n)横轴从1到32,纵轴从0到1.5 X=fft(x); %X(n)为x(n)的快速傅里叶变换

H=fft(h); %H(n)为h(n)的快速傅里叶变换

Y=X.*H; %Y(n)=X(n)*H(n)

%Y=conv(x,h);

y=ifft(Y); %y(n)为Y(n)的傅里叶反变换

subplot(3,1,3) %在第三个子窗绘图y(n)横轴从1到32,纵轴从0到6 stem(abs(y));

title('y(n=x(n)*h(n))');

axis([1,32,0,6]);

结果：

实验总结与思考

1、在较短的傅里叶变换中，FFT的计算速度与DFT相比不是很明显，序列计算长度越长，计算时间差距越大，FFT较快；

2、对于不同序列的较小长度的频谱分析可能会得到相同的频谱，适当加倍长度会避免这种情况的发生；

3、对同一序列的不同间隔的FFT变换，在满足奈奎斯特定律的情况下也会产生栅栏效应、频谱泄露、旁瓣效应等，采取适当的方法可以减弱这些不利效应；

4、在计算两个序列的离散卷积的时候要注意序列的长度L>=M+N-1。

实验三、IIR数字滤波器设计

源程序：

function pr7()

wp=0.2*pi; %频率转换（WP:通带截止频率）

ws=0.3*pi; %频率转换(WS:阻带截止频率)

Rp=1; %对于参数的说明

As=15;

T=1;

Fs=1/T;

OmegaP=(2/T)*tan(wp/2); %OmegaP(w)=2*tan(0.1*pi)

OmegaS=(2/T)*tan(ws/2); %OmegaS(w)=2*tan(0.15*pi)

ep=sqrt(10^(Rp/10)-1);

Ripple=sqrt(1/(1+ep.^2));

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

[cs,ds]=u_buttap(N,OmegaC);

[b,a]=bilinear(cs,ds,Fs);

[mag,db,pha,w]=freqz_m(b,a);

subplot(3,1,1); %统共要绘制三个图，把此图放在第一个子窗绘制

%幅度响应的图形

plot(w/pi,mag);

title('幅度响应')

xlabel('w(pi)');

ylabel('H');

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

set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);

set(gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1]);

grid;

subplot(3,1,2); %在第二个子窗以分贝为单位绘制幅度响应的图形plot(w/pi,db);

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

xlabel('w(pi)');

ylabel('H');

axis([0,1,-40,5]);

set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);

set(gca,'YTickmode','manual','YTick',[-50,-15,-1,0]);

grid;

subplot(3,1,3); %在第三个子窗绘制相位响应的图形plot(w/pi,pha);

title('相位响应');

xlabel('w(pi)');

ylabel('pi unit');

%axis([0,1,0,1.1]);

set(gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]);

set(gca,'YTickmode','manual','YTick',[-1,0,1]);

grid;

function [b,a]=u_buttap(N,OmegaC)

[z,p,k]=buttap(N);

p=p*OmegaC;

k=k*OmegaC.^N;

B=real(poly(z));

b0=k;

b=k*B;

a=real(poly(p));

function [mag,db,pha,w]=freqz_m(b,a)

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

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

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

mag=abs(H);

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

pha=angle(H);

结果：

实验总结与思考

1、数字滤波器可以用模拟滤波器通过模拟到数字的变换来实现，具体方法有三种：冲激响应不变法、阶跃响应不变法和双线性变换法等；

2、冲激响应不变法仅适合于基本上是限带的低通滤波器，对于高通或带阻滤波器应该附加限带要求，以避免严重的混叠失真；

3、双线性变换法可以用频率预畸变的方法来补偿频率标度的非线性失真；

4、数字低通滤波器的模拟截止频率是抽样时间的函数；

5、低通滤波器对相应频率的输入信号有相应限制作用。

实验四、FIR数字滤波器的设计

源程序：

function pr8()

wp=0.2*pi;

ws=0.3*pi;

tr_width=ws-wp;

M=ceil(6.6*pi/tr_width)+1;

n=0:M-1;

wc=(ws+wp)/2;

alpha=(M-1)/2;

m=n-alpha+eps;

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

w_ham=(hamming(M))';

h=hd.*w_ham;

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

delta_w=2*pi/1000;

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

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

subplot(2,2,1);

stem(n,hd);

title('理想冲激响应');

axis([0,M-1,-0.1,0.3]);

ylabel('hd(n)');

subplot(2,2,2);

stem(n,h);

title('实际冲激响应');

axis([0,M-1,-0.1,0.3]);

ylabel('h(n)');

subplot(2,2,3);

plot(w/pi,pha);

title('滤波器相位响应');

axis([0,1,-pi,pi]);

ylabel('pha');

set(gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]); set(gca,'YTickmode','manual','YTick',[-pi,0,pi]); grid;

subplot(2,2,4);

plot(w/pi,db);

title('滤波器幅度响应');

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

ylabel('H(db)');

set(gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]); set(gca,'YTickmode','manual','YTick',[-50,-15,0]);

function [mag,db,pha,w]=freqz_m(b,a)

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

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

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

mag=abs(H);

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

pha=angle(H);

结果：

实验总结与思考

1、IIR滤波器的系统函数是有理分式，其分母多项式对应于反馈支路，因而这种滤波器是递归结构的系统，只有当所有极点都在单位圆内时滤波器才是稳定的；

2、IIR滤波器可以用比FIR滤波器少的阶数来满足相同的技术指标，这样，IIR滤波器所用的存储单元和所需的运算次数都比FIR滤波器少；

3、IIR滤波器可利用模拟滤波器现成的设计公式、数据和表格，计算工作量较小，对计算工具要求不高；FIR滤波器没有现成的设计公式，窗函数法只给出窗函数的计算公式，但计算通带和阻带衰减仍无显式表达式，一般，FIR滤波器的设计只有计算机程序可以利用；

4、FIR滤波器的相位严格线性；

5、相同指标情况下，FIR滤波器的系统特性更好些；

6、改变FIR滤波器的阶数，对相位特性的影响比幅度特性更显著；

7、如果需要严格线性相位的滤波器，则需要用FIR滤波器。

综合实验

源程序：

function pushbutton1_Callback(hObject, eventdata, handles)

Fs=12000; %抽样频率

time=10; %录音时间

x=wavrecord(time*Fs,Fs,'double'); %录制语音信号wavwrite(x,Fs,' F:\52102426\20130829.wav '); %存储语音信号

function pushbutton2_Callback(hObject, eventdata, handles)

y=wavread(' F:\52102426\20130829.wav ');

wavplay(y,12000); %播放录音

function pushbutton3_Callback(hObject, eventdata, handles)

[z1,fs,bits]=wavread(' F:\52102426\20130829.wav ')

y1=z1(1:18000);

Y1=fft(y1);

wp=0.2*pi;

ws=0.25*pi;

wdel=ws-wp;

N=ceil(8*pi/wdel);

Wn=(ws+wp)/2;

window=hamming(N+1);

b=fir1(N,Wn/pi,window);

x=fftfilt(b,z1);

sound(x,fs,bits);

function pushbutton4_Callback(hObject, eventdata, handles)

[z1,fs,bits]=wavread(' F:\52102426\20130829.wav ')

y1=z1(1:18000);

Y1=fft(y1);

wp=0.2*pi;

ws=0.25*pi;

tr_width=ws-wp;

M=ceil(8*pi/tr_width)+1;%求序列长度

n=0:M-1;

wc=(ws+wp)/2;

alpha=(M-1)/2;

m=n-alpha+eps;%右移

hd=sin(wc*m)./(pi*m);%给定的HD

w_ham=(hamming(M))';%汉明窗

h=hd.*w_ham;%加窗

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

axes(handles.axes1);

plot(w/pi,db);

title('滤波器幅度响应');

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

ylabel('H(db)');

function [mag,db,pha,w]=freqz_m(b,a)

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

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

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

mag=abs(H);

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

pha=angle(H);

function pushbutton5_Callback(hObject, eventdata, handles)

[z1,fs,bits]=wavread(' F:\52102426\20130829.wav ')

y1=z1(1:18000);

Y1=fft(y1);

wp=0.2*pi;

ws=0.25*pi;

wdel=ws-wp;

N=ceil(8*pi/wdel);

Wn=(ws+wp)/2;

window=hamming(N+1);

b=fir1(N,Wn/pi,window);

x=fftfilt(b,z1);

X=fft(x,18000);

axes(handles.axes4);

plot(abs(Y1));axis([0,10000,0,100]);

title('滤波前信号频谱');

axes(handles.axes2);

plot(z1);

title('滤波前信号波形');

axes(handles.axes5);

plot(abs(X));axis([0,10000,0,100]);

title('滤波后信号频谱')

axes(handles.axes3);

plot(x);

title('滤波后信号波形');

实验结果及分析

1、录音抽样频率设定为12000Hz，录音时间设置为5s，滤波器采用FIR数字低通滤波器，使用汉明窗设计，通带截止频率为0.1kHz，阻带截止频率为0.15kHz，根据要求设计出的滤波器满足指标要求。

2、按钮1为录音键、按钮2为播放录音文件按键、按钮3为LPF滤波后播放按键、同时画出图形滤波器幅度响应，相对应上面的为滤波器幅度响应曲线输出窗口、滤波前后所录声音文件的波形和频谱。