数字信号处理三级项目

1、数字信号处理三级项目14-通信工程4班目录目录01时域信号的采样及谱分析Part One02设计FIR滤波器Part Two03利用FIR滤波器滤波Part Three04思考题Part FourCONTENTS01时域信号的采样及谱分析Part One谱分析过程2.对截取的模拟信号进行采样（选择满足时域采样定理的频谱采样率）1.对模拟信号进行截取（选择合适的窗函数）3.对采样所得序列进行DFT(FFT)运算4324.信号的频谱分布特性1利用绘制时域采样序列图及频谱图利用+easyx绘制时域采样序列图及频谱图利用绘制时域采样序列图及频谱图信号的频谱分布特性取100点，补零至128点取100点，

2、补零至256点，256点DFT高密度谱取256点，256点DFT高分辨率谱取100点，补零至256点，256点DFT运算量 一般来说,FFT比DFT运算量小得多,N点的FFT需要做(N/2)log2N次乘法运算,而N点DFT需要做N2次乘法运算,由此看来N点DFT运算量大约是FFT的2N/log2N倍,例如对1 024点的变换,DFT大约是FFT的200倍.然而实际应用时存在下列情况: 点数或采样率的可选性 对DFT来讲,其变换点数可任意选定,.FFT的变换点数必须是有规律的 实时性 DFT运算可以用采一点后立即进行相乘、累加运算的方法FFT运算必须在全部点采集结束后才能开始进行计算,因此从某

3、种角度讲DFT的实时性优于FFT.DFT与FFTCONTENTS02设计FIR滤波器Part Four窗函数法频率取样法切比雪夫等波纹逼近法滤波FIR滤波器的主要设计方法01020304利用窗函数设计FIR滤波器根据对阻带衰减及过渡带的指标要求，选择窗函数的类型，并估计窗口长度N先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣窄的窗函数。然后根据过渡带宽度估计窗口长度N。构造需要逼近的频率响应函数计算若已知通带边界频率和阻带边界频率，则取加窗得到设计结果：)( nhdeHeHNjwdgjwdw2/ )1()()()()()(nwnnhhd2wwwspcdwnee

4、Hhjwnjwdd)(21)(常用的窗函数的幅频特性矩形窗三角形窗汉宁窗海明窗布莱克曼窗窗函数设计FIR低通滤波器的频率特性矩形窗矩形窗三角形窗三角形窗汉宁窗汉宁窗海明窗海明窗布莱克曼窗布莱克曼窗不同窗函数对滤波性能的影响窗函数类型窗函数类型旁瓣峰值旁瓣峰值n/dB过渡带宽度过渡带宽度Bt阻带最小衰减阻带最小衰减s/dB近似值精确值矩形窗-134/N1.8/N-21三角窗-258/N6.1/N-25汉宁窗-318/N6.2/N-44海明窗-418/N6.6/N-53布莱克曼窗-5712/N11/N-74窗函数的选择标准1、较低的旁瓣幅度2、旁瓣幅度要下降得快，以利于增加阻带衰减。3、主瓣宽度要

5、窄，这样滤波器过渡带较窄。CONTENTS03利用FIR滤波器滤波Part five之之3.1.滤波前后的频域波形（采样点数为200点）各滤波器截止频率选择低通滤低通滤波器波器高通滤高通滤波器波器带通滤带通滤波器波器带阻滤带阻滤波器波器fc(Hz)204211fc1=204fc1=204fc2=211fc2=211c(rad)0.4080.4220.4080.408矩形窗低通&高通矩形窗矩形窗低通低通FIR滤波器滤波器矩形窗矩形窗高通高通FIR滤波器滤波器矩形窗带阻&带通矩形窗带阻FIR滤波器矩形窗带通FIR滤波器其他窗函数的低通滤波器汉宁窗汉宁窗低通低通FIR滤波器滤波器海明窗海明窗低通低通

6、FIR滤波器滤波器其他窗函数的高通滤波器汉宁窗汉宁窗高通高通FIR滤波器滤波器海明窗海明窗高通高通FIR滤波器滤波器其他窗函数的带阻滤波器汉宁窗汉宁窗带阻带阻FIR滤波器滤波器海明窗海明窗带阻带阻FIR滤波器滤波器其他窗函数的带通滤波器汉宁窗汉宁窗带通带通FIR滤波器滤波器海明窗海明窗带通带通FIR滤波器滤波器之之3.2.提高采样时间，对信号滤波前后时域频域分析滤波前的时域波形取采样点为取采样点为1000点滤波点滤波前的时域图前的时域图低通滤波器滤波后的滤波后的时域图时域图低通滤低通滤波器波器fc(Hz)204c(rad)0.408海明窗的设计，阶数，阶数N=475带通滤波器滤波后的滤波后的时

7、域图时域图带通滤带通滤波器波器fc(Hz)fc1=204fc2=211c(rad)0.408海明窗的设计，阶数，阶数N=475带阻滤波器滤波后的滤波后的时域图时域图带阻滤带阻滤波器波器fc(Hz)fc1=204fc2=211c(rad)0.408海明窗的设计，阶数，阶数N=475高通滤波器滤波后的滤波后的时域图时域图高通滤高通滤波器波器fc(Hz)211c(rad)0.422海明窗的设计，阶数，阶数N=475CONTENTS04思考题Part Two用DFT对信号进行频谱分析时，f1=201hz,f2=208hz,f3=214hz,采样频率fs=1000hz,观测时间为Tp=0.1s时的时域和

8、频域波形：频域波形 时域波形采样点N=100，补0至128此时采样分辨率为F=1/Tp=1/0.1s=10hz,而f1、f2、f3之间互相差7hz、6hz，10hz不是它们的整数倍，所以不能将三个单倍频分辨出来。保持信号长度M不变，DFT变换的N=256时信号波形：频域波形时域波形采样点N=256此时相对于N=128时频谱更加密集,其包络更加接近序列的FT。在分辨率F=1/Tp=1/0.1s=10hz,则在保持观测时间不变的条件下不能提高分辨率。此时由f/fs=k/N可算得k1=20.1,k2=20.8,k3=21.4采样点为20点（整周期）采样点为25点（非整周期）采样点为30点（非整周期）

9、采样点为80点（整周期）当采样点N=32,并进行32点DFT，经过计算k=5,27的谱线代表f=40Hz； k=10,22的谱线代表f=80Hz；但对于f=60Hz，对应的k=7.5,60Hz以频谱泄露的形式出现80Hz当采样点N=64,并进行64点DFT，经过计算k=10,54的谱线代表f=40Hz； k=20,44的谱线代表f=80Hz； k=15,49的谱线代表f=60Hz；40Hz80Hz60Hz40Hz当输入频率f=k(Fs/N)，k为整数，相当于时域被周期截取，第k根谱线就表示该频率；当k为小数时，相当于时域被非周期截取，发生频谱泄露，所有谱线都对该频率有响应。频域波形时域波形采样

10、频率fs=500hz时的波形图当采样频率降为500hz时，频率分辨率F=fs/N=500/100hz=10hz，此时仍旧不能将三个单频信号分辨出来。频域波形时域波形Tp=0.2s时的波形当观测时间为0.2s时，F=1/Tp=1/0.2s=5hz,f2、f3之间互相差7hz、6hz，5hz不是它们的整数倍,此时不能将三个单频信号分辨出来。 改变输入信号各个频率分量的频率，f1=199hz,f2=200hz,f3=204hz,fs=1000hz,Tp=0.2s时的输出波形：频域波形时域波形f1=201hz,f2=208hz,f3=214hz,采样频率fs=1000hz,观测时间为Tp=0.2s时输

11、出波形：对比两次频谱输出，可知，改变输入信号频率分量之后，输出的频域波形在中心频率附近的旁瓣幅度增强。若要正确对该信号进行频谱分析，则需要提高观测时间Tp。利用快速卷积对信号进行滤波需注意：1、要求对短序列补很多零点，而且长序列必须全部输入后才能进行快速卷积计算。2、在采样前进行预滤波，滤出高于折叠频率FS/2的频率成分，以免发生频谱混叠现象3、由于栅栏效应，有可能漏掉大的滤波分量。为了把原来被漏掉的频谱分量检测出来，对有限长序列，可以再原序列尾部补零；无限长序列，可以增大截取长度及DFT变换区间长度。4、实际遇到的序列可能是无限长的，用DFT对其进行谱分析时，必须将其截短，形成有限长序列称为

12、窗函数。FIR数字滤波器具有线性相位的充要条件是偶对称 即h(n)=h(N-1-n)或奇对称 即h(n)=-h(N-1-n)滤波器具有线性相位的条件采样点为20点（整周期）采样点为25点（非整周期）采样点为80点（整周期）波形两端突变变得平滑些以此压低谱窗的旁瓣。因为旁瓣泄露量最大旁瓣小了泄露也相应减少了。加矩形窗对滤波器频率响应的影响参考文献 1 高西全、丁玉美.数字信号处理（第三版）M.西安：西安电子科技大学出版社，2015 2 高西全、丁玉美.数字信号处理（第三版）M学习指导.西安：西安电子科技大学出版社，2015 3 王大伦.数字信号处理.北京：清华大学出版社M，2014 4 Sanj

13、it K.Mitra.数字信号处理基于计算机的方法（第三版）M.北京：电子工业出版社，2006附录程序代码#include #include#includedrawfunction.h#includeDFS.husing namespace std;complex aa3000=0;float h10000=0;float A110000=0;float B110000=0;float A210000=0;float B210000=0;float a10000=0;float b10000=0;float IMM10000=0;float REE10000=0;float C10000=0;

14、float LC10000=0;complex bb3000;complex cc3000;int N1=200;/*选择滤波器的种类*/void selectwf()int k;printf(Please select the Wave filter:n);printf(1-Low Pass Wave Filer(低通滤波器)n);printf(2-High Pass Wave Filer(高通滤波器)n);printf(3-Stop Band Wave Filer(带阻滤波器)n);printf(4-Pass Band Wave Filer(带通滤波器)n);scanf(%d,&k); s

15、witch(k) case 1: DT(b,N);break; case 2: GT(b,N);break; case 3: DaiZ(b,N);break; case 4: DaiT(b,N);break; printf(Please select the Window Function:n); printf(1-RECTANGLE WINDOW(矩形窗)n); printf(2-HANNING WINDOW(汉宁窗)n); printf(3-HANMING WINDOW(哈明窗)n); printf(4-BLACKMAM WINDOW(布莱克曼窗)n); scanf(%d,&k); swi

16、tch(k) case 1: JXC(B1,N);break; case 2: HNC(B1,N);break; case 3: HMC(B1,N);break; case 4: BLKMC(B1,N);break; for(int i=0;iN;i+) hi=bi*B1i; /*进行滤波*/void startfilter() Convolution(a,h,C,N1,N); for(int i=0,j=(N-1)/2;iN1;i+) Ci=Cj+i; DFT(C,A1,A2,REE,IMM,N1);/*图形的绘制*/void graph() DFT(a,A1,LC,REE,IMM,N1);

17、 initgraph(1400, 1000, SHOWCONSOLE | NOCLOSE); show(LC,N1,50,1400,300,150,450,50); /show(C,N1,50,1400,300,150,450,50); show(a,N1,50,1400,500,400,600,50); /show(A2,N1,50,1400,500,400,600,50); /show(C,N1,50,1400,700,550,700,50); /show(a,N1,50,1400,700,550,700,50); getchar(); system(pause);void main( )

18、 quyang(a,1000,N1); selectwf(); startfilter(); graph();#includecomplex/#includedrawfunction.husing namespace std;#define pi 3.1415926double wc=0.408*pi;/204HZ=0.408，211HZ=0.422int N=475;void Convolution(float *Tab_A,float *Tab_B,float *cout,int la,int lb) int ii=0,jj=0; int k=0;float out10000=0;for(

19、 int p=0;p=la+lb-2;p+) /卷积后的序列号 for(int i=0;i=la-1;i+) /序列1的序号 for( int j=0;j=lb-1;j+) /序列2的序号 if(i+j)=p)outk=outk+Tab_Ai+abs(ii)*Tab_Bj+abs(jj); /printf(%dn,outk); /输出卷积后的序列 coutk=outk; k=k+1; /*DFT函数*/*a为输入序列x(n)输入可为复数序列，c记录得到的|X(K)|,c输入必须为0*/*re记录X(K)的实部，im记录X(K)的虚部，N为周期的大小*/void DFT(float *a,flo

20、at *b, float *c,float *re,float *im, int N) complex po10000; float op10000;int n;float s = float(-(2 * pi) /(float)N);complex tm; for(n=0;nN;n+) pon=complex(an, bn);/使得输入序列为复数 for(n=0;nN;n+)opn=1;for (int k = 0; k N; k+)tm = complex(ck, 0);for ( n = 0; n N; n+)tm += pon*polar(opn,s*k*n); rek=real(tm

21、); imk=imag(tm); ck =abs(tm);/printf(%ft,ck);/*取样函数*/void quyang(float *Out,float fs,int length) float f1=201; float f2=208;float f3=214;/*float f1=201; float f2=208;float f3=214;*/ for(int i=0;ilength;i+)/Outi=(float)cos(2*pi*f1/fs*i)+(float)cos(2*pi*f2/fs*i)+(float)cos(2*pi*f3/fs*i);Outi=(float)co

22、s(2*pi*f1/fs*i)+(float)cos(2*pi*f2/fs*i)+(float)cos(2*pi*f3/fs*i);/*各类窗函数*/*矩形窗*/void JXC(float *wn,int n)int i;for(i=0;in;i+) if(iN)wni=1;elsewni=0;/*三角窗*/void SJC(float *wn,int n)int i;for(i=0;i(N-1);i+)if(i=1/2*(N-1) wni=2*n/(N-1);elsewni=2-2*n/(N-1);for(i=N-1;in;i+)wni=0;/*汉宁窗*/void HNC(float *w

23、n,int n)int i;for(i=0;in;i+) if(iN)wni=(0.51-0.46*cos(2*pi*i/(N-1);elsewni=0;/*哈明窗*/void HMC(float *wn,int n)int i;for(i=0;in;i+) if(iN) wni=0.5*(1-cos(2*pi*i/(N-1);elsewni=0;/*布莱克曼窗*/void BLKMC(float *wn,int n)int i;for(i=0;in;i+) if(iN) wni=0.42-0.5*cos(2*pi*i/(N-1)+0.08*cos(4*pi*i/(N-1);elsewni=0

24、;/*各类理想滤波器*/*理想低通滤波器*/void DT(float *Hd,int n) int i; float a; a=(N-1)/2; for(i=0;iN;i+) if(i=a)Hdi=wc/pi; else Hdi=sin(wc*(i-a)/float(pi*(i-a); for(i=N;in;i+) Hdi=0; /*理想高通滤波器*/void GT(float *Hd,int n)int i;int a;a=(N-1)/2; for(i=0;iN;i+) if(i=a)Hdi=1-0.422; else Hdi=sin(pi*(i-a)/(pi*(i-a)-sin(0.42

25、2*pi*(i-a)/float(pi*(i-a); for(i=N;in;i+) Hdi=0; /*理想带通滤波器*/void DaiT(float *Hd,int n)int i;int a;a=(N-1)/2;for(i=0;iN;i+) if(i=a)Hdi=-0.408+0.422;else Hdi=-sin(0.408*pi*(i-a)/float(pi*(i-a)+sin(0.422*pi*(i-a)/float(pi*(i-a); for(i=N;in;i+) Hdi=0;/*带阻滤波器*/void DaiZ(float *Hd,int n)int i;int a;a=(N-1

26、)/2;for(i=0;iN;i+) if(i=a)Hdi=1-0.422-0.408; /else Hdi=sin(pi*(i-a)/(pi*(i-a)-sin(0.408*pi*(i-a)/float(pi*(i-a)-sin(0.422*pi*(i-a)/float(pi*(i-a);else Hdi=sin(pi*(i-a)/(pi*(i-a)-sin(0.422*pi*(i-a)/float(pi*(i-a)-sin(0.408*pi*(i-a)/float(pi*(i-a); for(i=N;in;i+) Hdi=0;#includegraphics.h/*绘图函数*/*输入为：序

27、列、序列长度、X轴的起始x像素点，X轴的末尾x像素点、X轴的y像素点*/*Y轴的x像素点、Y轴的起始y像素点、Y轴的末尾y像素点*/*绘图界面初始化函数需要在主函数中编写*/void show(float *AL,int Length,int X_x1=200,int X_x2=600,int X_y=200,int Y_y1=100,int Y_y2=300,int Y_x=200) /*TCHAR rs1 = _T(输入序列1为R5=1,1,1,1,1); outtextxy(500, 275, rs1); TCHAR rs2 = _T(输入序列2为R10=1,1,1,1,1,1,1,1,

28、1,1); outtextxy(500, 300, rs2);TCHAR ss = _T( ); outtextxy(50, 175, ss); TCHAR ss1 = _T( ); outtextxy(50, 475, ss1);*/initgraph(1000, 1000, SHOWCONSOLE | NOCLOSE);/int L;/L=Length; int i=0;int j=0; float dx=5;float dy=5;float djx=10;float djy=10;float ddx=0;float ddy=0;float max=0;line(Y_x,Y_y1,Y_x,

29、Y_y2);/Y轴line(Y_x-djx,Y_y1+djy,Y_x,Y_y1); /箭头line(Y_x+djx,Y_y1+djy,Y_x,Y_y1);line(X_x1,X_y,X_x2,X_y);/X轴line(X_x2-djx,X_y-djy,X_x2,X_y);/箭头line(X_x2-djx,X_y+djy,X_x2,X_y);ddx=(X_x2-Y_x)/Length; TCHAR s3; _stprintf(s, _T(%d), Length); / 高版本 VC 推荐使用 _stprintf_s 函数 / outtextxy(X_x1+(X_x2-X_x1)/3+130,Y_

30、y1-20, s); outtextxy( X_x1+2*(X_x2-X_x1)/3-10,Y_y1-3, s); /*设置X轴的坐标和标度*/for(i=0;i500) if(i%50=0)_stprintf(s, _T(%d), i); outtextxy(X_x1+i*ddx+2,X_y+2, s);else else if(Length100) if(i%20=0)_stprintf(s, _T(%d), i); outtextxy(X_x1+i*ddx+2,X_y+2, s);else else if(Length10) if(i%5=0) _stprintf(s, _T(%d),

31、i); outtextxy(X_x1+i*ddx+2,X_y+2, s); else else _stprintf(s, _T(%d), i); outtextxy(X_x1+i*ddx+2,X_y+2, s); /*找出序列中最大的值*/for(i=0;imax)max=ALi;float H=max;max=max*9/7;ddy=(X_y-Y_y1)/max;/*设置Y轴的坐标和标度*/for(i=1;i100)if(i%10=0) line(Y_x,X_y-i*ddy,Y_x+dx,X_y-i*ddy);/y坐标设置 line(Y_x,X_y+i*ddy,Y_x+dx,X_y+i*dd

32、y);/y坐标负轴 else line(Y_x,X_y-i*ddy,Y_x+dx,X_y-i*ddy);/y坐标设置 line(Y_x,X_y+i*ddy,Y_x+dx,X_y+i*ddy);/y坐标负轴if(H1000)if(i%50=0) line(Y_x,X_y-i*ddy,Y_x+dx,X_y-i*ddy);/y坐标设置if(H1000)if(i%100=0) _stprintf(s, _T(%d), i); / 高版本 VC 推荐使用 _stprintf_s 函数 outtextxy(Y_x-35,X_y-i*ddy, s); else else if(H100) if(i%20=0

33、) _stprintf(s, _T(%d), i); / 高版本 VC 推荐使用 _stprintf_s 函数 outtextxy(Y_x-30,X_y-i*ddy, s); else else if(H10) if(i%5=0) _stprintf(s, _T(%d), i); / 高版本 VC 推荐使用 _stprintf_s 函数 outtextxy(Y_x-20,X_y-i*ddy, s); _stprintf(s, _T(%d), -i); / 高版本 VC 推荐使用 _stprintf_s 函数 outtextxy(Y_x-25,X_y+i*ddy, s); else else _

34、stprintf(s, _T(%d), i); / 高版本 VC 推荐使用 _stprintf_s 函数 outtextxy(Y_x-10,X_y-i*ddy, s); _stprintf(s, _T(%d), -i); / 高版本 VC 推荐使用 _stprintf_s 函数 outtextxy(Y_x-15,X_y+i*ddy, s); /*y的负半轴*/*绘制序列点*/ setfillcolor(WHITE);for(i=0;iLength;i+) fillcircle(X_x1+i*ddx,X_y-ALi*ddy,2); line(X_x1+i*ddx,X_y,X_x1+i*ddx,X

35、_y-ALi*ddy);/*修改其中的汉字即可显示不同的值*/ setbkcolor(GREEN);/ cleardevice(); setcolor(WHITE);RECT r1=X_x1+(X_x2-X_x1)/3,Y_y1-25,X_x1+2*(X_x2-X_x1)/3,Y_y1;drawtext(_T(序列点数为：), &r1, DT_CENTER | DT_VCENTER | DT_SINGLELINE); TCHAR s1 = _T(通信4班2组：); outtextxy(650, 10, s1); TCHAR s2 = _T(组长：张曌); outtextxy(660, 30, s2); TCHAR s3 = _T(组员：郑岚方); outtextxy(660,50, s3); TCHAR s4 = _T(组员：蒲琼); out