数字信号处理实验

1DFTDFTDFTDFT(02的Nk2kN(0kN1上的N2kX(ej)|2k

N k

j2kn

X 0kN (2x(n的NDFTX(kx(n在Nk2kN(0

N1X(k1：X(kX(ej2.12.1.NDFTDTFT2.1X(ej)

x(n)e

N

X(k

kne

(22-2

nk NN

N其中(x

2

方法2：然而在实际计算中，上诉插 不见得是最好的方法DFTDTFT2

果我们增加数据的长度N，DFTDTFT的结果，这样可以利用DFT来近似计算DTFT。如果没 的数据，可xt，按采样间隔TM，X(j)

Mx(t)ejtdt

(nT)e

(2

X(j进行NX(j)

M2

j2kn

TXM

DFT确定时域采样间隔T，x(n确定截取长度M，得到MxM(n)x(n)w(n)w(n确定频域采样点数N，NMFFTNDFT，XM(k)根据式（2-6）XM(kX(j采样点的近似值。DFTN实验中的DFT运算可以采用中提供的FFT来实现。DTFT可以利用矩阵运算的方法进行计算。x(n)21,1,1}DTFT，并画出[,]4DFT，并把结果显示在（1）x(n64DFT，DFTDTFT，DTFT，并画出[,]n=0:3;x=[2-111];w=->>axistightxlabel('w');title('Phase');%画出相位谱axistight321 0

w4DFT，并把结果显示在（1）n=0:3;x=[2-11>>subplot(211);>>holdstem([0:pi/2:3*pi/2],abs(Y));4>>>>axis>>subplot(212);plot(w,angle(X)/pi);holdon>>stem([0:pi/2:3*pi/2],angle(Y4>>xlabel('w');title('Phase');axistight3210 w10 wx(n64DFT，并显示结果。n=0:3;x=[2-111];Y=fft(x,64);%64subplot(211);axistightaxistight3210

10 kDFTDTFT，（2（3）DFTDTFT序列后面补零的方法，则可以提高采样密度，使DFT的点数，从而得到精DTFTx(n)cos(0.48n)cos(0.520n10DFTx(nx(n100DFTx(n0n100DFTx(n0n10DFTx(nx(n100DFTx(n)的频谱。要求画出相应波形。0n10DFTx(n>>>>subplot(211);>>xlabel('0<=n<=10,DFTx(n)>>axis>>xlabel('0<=n<=10,DFTx(n)频谱');title('Phase');axistight8642 0<=n<=10,DFT估计x(n)频谱10 0<=n<=10,DFT估计x(n)频谱x(n100DFTx(n>>>>subplot(211);DFT>>axis>>xlabel('0<=n<=10,100，DFTx(n)频谱');title('Phase');axistight8642 0<=n<=10,补零至100点，DFT估计x(n)频谱10 0<=n<=10,补零至100点，DFT估计x(n)频谱0n100DFTx(n>>>>subplot(211);>>xlabel('0<=n<=100，DFTx(n)>>axis>>subplot(212);xlabel('0<=n<=100，DFTx(n)频谱');title('Phase');axistight 0<=n<=100，DFT估计x(n)频谱210 0<=n<=100，DFT估计x(n)频谱x(t0.15sin(2f1tsin(2f2t0.1sin(2f3t)f11Hzf22Hzf33Hzx(t的表达式可以看出，它包含三个频率的正弦波，但DFTx=0.15*sin(2*pi*t)+sin(2*pi*2*t)->>>>subplot(211);>>xlabel('DFTx(n)频谱>>axis>>xlabel('DFTx(n)频谱');title('Phase');axistight0 DFT估计x(n)频谱20 DFT估计x(n)频谱DFT（幅度谱。分析采用不同的采样间0120，266543210 k当然在实验中也遇到不少，由于编程能力有限，我对一些函数的用法不甚了解，在2FFT掌握利用FFT加深理解相加法和保留法利用FFT设x(n)为L点序列，h(n)为M点序列，x(n)h(n)y1(n)x(n)*h(n)

x(m)h(n

（3-y1(n的长度为L+M-1x(n)和h(n)y(n)

Nh(n) x(m)h(nm)NRN

（3-(3-DFT性质，x(n)h(n)的NDFTDFT(3-4)FFT算法必须选择；为了能使用基-2算法，要求 法使得x(n)和h(n)的长度均为N。计算x(n)和h(n)的N点FFTx(n)

X(k

h(n)

HY(k)X(k)H利用IFFTY(k)IDFTY(k)

我们单位取样响应为h(n)的线性系统，输入为x(n)，输出为y(n)，以细分为保留法和相加法。保留法：设x(n)的长度为，h(n)的长度为M。把序列x(n)分成多段点序列，每段雨前一段重写M-1个样本。并在第一个输入段前面补M-1去，只保留后面N-M+1个正确的输出样本，把它们合起来得到总的输出。利用FFT实现保留法的步骤如下在x(n)前面填充M-1将x(n)分为若干段N点子段，设L=N-M+1为每一段的有效长度，则第i段的计算每一段与h(n)的N点圆周卷积，利用FFTx(n) X(k h(n)

HY(k) X(k)H iY(k)i

舍去每一段卷积结果的前M-1个样本，连接剩下的样本得到卷积结果y(n)。x(n)xi每一段卷积的长度为L+M-1，所以在做求和时，相邻两段序列由M-1个本，即前一段的最后M-1个样本和下一段前M-1个样本序列，这个重叠部分相加，再与不的部分共同组成y(n)。利用FFT实现保留法的步骤如下将x(n)分为若干L点子段计算每一段与h(n)的卷积，根据快速卷积法利用FFT将各段相加，得到输出y(n)y(n)yi(n4、可能得到的函实验中FFT运算可采用中提供的函数fft来实现x(n)=u(n)-u(n-L),0≤n≤Lh(n)=cos(0.2πn),0≤n≤ML=M,根据线性卷积的表达式和快速卷积的原理分别编程实现计算两个序列L=4096M=256编程实现利用相加法计算两个序列的线性卷积,L=2048且M=256时2编程实现利用保留法计算两个序列的线性卷积,L=2048且M=256时2fori=1:L-x=[1h=[h时间已过0.004177时间已过0.002205时间已过0.008691时间已过0.014703时间已过0.021067时间已过0.036894时间已过0.058499时间已过0.080003结论：在LMfori=1:L-x=[1h=[h时间已过0.0078410.009291时间已过0.050958时间已过0.080958forj=1:L-x=[1h=[hx=[x,zeros(1,N-t=zeros(1,M-1);H=zeros(1,L+M-1);a=floor(L/N);fork=0:aH1=fft(A,L+M-1);H2=fft(h,L+M-1);HH=H1.*H2;t(1:M-1)=q(N+1:N+M-0.003830fori=1:L-x=[1y=[ya=floor((L+M-2)/(LL))+1;fork=0:a-1b=fft(xk,N);C=fft(h,N);Z=b.*C;0.007472，DFTDFT.3IIR器的单位冲激响应ha(t，h(n)等于ha(t的取样值。率响应，进而求得数字滤波器的系统函数H（zH(s变换时域采样h(nT

H 将Ha(sHa(ssk N

(s进行拉式变换h(t)Aepk k 对h(th(n)h(n)Aepkntu(nT)A(epkT)naNh(n)zH(zN

k

k

k11epkT通过正切变换Ω=tan(Ω1T/2)sjΩs1jΩ1jπ/Tjπ/Tc=2/T。确定数字滤波器的性能指标p,stRp,Asp

根据指标pstRpAsHa(s Ha(sHa(s)=spk kp转换成数字极点epkTk kH(z)1epTz1kkHa(sH(z)1s

1eskT

z1：residueresiduezs

sMb sM1 bH(s)M M1 sN sN1

H(s)

s

s

bbz1 bzMH(z) aas1 asNH(z)

1p1pkk(z) N1p 1p 2：（1）中提供了impinvar函数采用脉冲响应不变法实现模拟滤波器到向量babzaz，fs（1[bz,az]=impinvar(b,a1（2）中提供了bilinear函数采用双线性变换法实现模拟滤波器到数babzaz，fs设采样fs=10kHz，设计数字低通滤波器，满足如下指标阻带截止频率：fst=1.5kHz，阻带衰减IIR[NOmegaC]=buttord(Wp,Ws,Rp,As,'s');[ba]=butter(N,OmegaC,'s');[bzaz]=impinvar(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');10

1PhasePhaseof0

50 切比Ⅰ型：[N[bzaz]=impinvar(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');0

1PhasePhaseof0

GroupGroup50 切Ⅱ型：[N[bzaz]=impinvar(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');0

1PhasePhaseof0

50 [N[bzaz]=impinvar(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');0

1PhasePhaseof0

0 [NOmegaC]=buttord(Wp,Ws,Rp,As,'s');[bzaz]=bilinear(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');0

1PhasePhaseof0

GroupGroup50 切Ⅰ型：[N[bzaz]=bilinear(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');10

1PhasePhaseof0

GroupGroup50 切Ⅱ型：[N[bzaz]=bilinear(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');0

1PhasePhaseof0

Group5Group0 [N[bzaz]=bilinear(b,a,10000);[Hw]=freqz(bz,az);subplot(221);plot(w/pi,abs(H));gridon;subplot(222);plot(w/pi,20*log10(abs(H)/max(abs(H))));gridon;subplot(223);plot(w/pi,angle(H)/pi);gridon;xlabel('\omega(\pi)');ylabel('PhaseofH(e^j^\omega)(\pi)');subplot(224);plot(w/pi,grd);gridon;xlabel('\omega(\pi)');ylabel('Groupdelay');10

1PhasePhaseof0

GroupGroup50 01kHz1dB；1.5kHz15dB。切I型低通滤波器满足题目要求在0到1kHz时满足通带波动1dB；1.5kHz15dB。切II型低通滤波器不满足题目要求，通带波动明显大于1dB，阻带衰15dB。15dB。01kHz1dB；1.5kHz15dB。切I型低通滤波器满足题目要求在0到1kHz时满足通带波动1dB；1.5kHz15dB。切II型低通滤波器不满足题目要求椭圆型低通滤波器不满足题目要IIR滤波器的问题。在斯低通滤波器、切I型低通滤波器、切II型低通滤波器、椭圆低通滤波器的不同之处理论知识的学实践中得到巩固在分析实验中的问题时DSP的学习不能生搬硬套，而应该自己多动手，多时4FIRFIRdH(ejd然后用窗函数截取它的单位脉冲响应hd(n)FIRdH(ejd相位，在阻带上具有零响应。一个带宽为c(c的低通滤波器由下式给H(ej)

c其中sinc(nhd(n)

(n

为了得到一个h(n长度为NFIRN 2用窗函数截取hd(n)得到所设计FIR数字滤波器：h(n)hd （Hamming）窗、布莱克曼（Blackman）窗、凯瑟（Kaiser）窗等。提7-1所示。表7-1中产生窗函数函数函数N数(n)7-27-2窗I1(12n)20 N (n) I0[I0x是修正的零阶贝塞尔函数，参数于相同的N可以提供不同的过渡带宽。由于贝塞尔函数比较复杂，这种窗函数已知给定的指标p,stRp和

，滤波器长度NstNAs7.9510.5842(A21)040.07886(A21),21

H(ejwH(k）,将Hd(k）H（k H(k）=H H(k）

h(n)=IDFT[H（k）],H(ejw)=FTRH（k)的幅度和相位一定满足线H(k)表示成如下形式：rrh(n)H(k)=H*(N-k),由此得到：Hr(k)=Hr(N-1FIRp0.2,Rpst0.3,Asw_boxcar=boxcar(N)';%hd(n)h(n)axistight;axistight;axistight;grid10 n

0 0 0 n

w_hanning=hanning(N)';hd(n)h(n)axistight;axistight;axistight;grid10 n

0 0

n

w_hamming=h