C语言实现数字信号处理算法.doc

中I屮:宁技术丿、宁电JT.F旳信息科学广媒体通伫3、 (Copyright 2000)附录A C语言实现数字信号处理算法附录A1 BC下复数类型的实现木部分内容可以在 hg:/202.3X.75.33/ds)/C/Com【lcxUsc.cpp 获得。1、利用BC提供的复数支持/BC中便用复数类型使用小例(ComplexUse.Cpp文件、林include ciosireain.h林include <complcx.h>int main(void)double x = 3.1. y = 4.2;complex z = complex(x.y):cout« ”z = u« z « W;cout « ” and imaginary real part = ” « imag(z)« "W；cout « "z has complex conjugate = M « conj(z) « * n”；return 0;2、定义复数类，填写相应成员函数/C中的芟数类型调用时可能不是年帘好用,叮门U怎义复数熬ComphxUse.Cpp文件)class Complex!public:Complex()Coinplex( float re. float im );float r()return real:;float i()return imag;float inod() return sq rt (re al *real+i mag *i n)ag);Complex operator+( Complex &other ):Complex operator Complex &othcr);Complex operator4c( Complex &other );Complex operator/( Complex &olher );private:float real, imag: J/ Operator overloaded using a member fitnetionComplex:Complex( float re.float iin)real=re:imag=im:;Complex Complex:operator+( Complex &other )return Complex( real + other.real, imag + other.imag );Complex Complex:operator Complex &other)return Complex( real - other.real. imag other.imag );Complex Complex:operator*( Complex &other )float x.y;x=real*other.real-imag*other.imag;y=real other.i mag+imag *other. real;return Complex( x,y):Complex Complex:operator/( Complex &olhcr )float x.yj;l=othcr.realother.real+otl)er.iinag*other.imag:x=real*other.real+imag*other.imag:y=other. real*i mag- real *other.i mag:x=x/l;y=y/l;return Complex(x.y);附录A2 BC下的绘图木部分内容可以右 htp:/2O2.38.75.33/dsp/C/PStBYBC.cpp 获得。1、通用绘图函数的有关说明(1) 函数输入参数说明left :绘图区的左上角横坐标top :绘图区的左上角纵坐标right:绘图区的右卜角横坐标bottom :绘图区的右卜角纵坐标f:需要绘制图形的数组length: f数组的长度(2) 函数调用示例在以(10.5)为庄上角坐标，(200.240)为右下角坐标的区域内绘制数组x的图形，数组氏度为10Plot(10.5,200,240.xJ0);BC图形程序盂耍“BC安装hl录bgiegavga.bgi"文件支持如果您在阅读过程M1发现疏涓和iJi误.诘您尽快和编者取得联系nclwork cxh中国科学技术人宁电f息科学 广媒体通伫®、 (Copyright 2000)2、通用绘图函数void Plot(int leftjnt lop,ini right.int bottontdouble %int length)int n:int Left.Top.RightBottom.Mid:double x_unil_lcngih.y_unil_lcngih.max= 1.0;for(n=begin:n<=end;n+) /begin和end我示数纟H F标的起始值和结束值 if(max<fabs(Hn|)max=fabs(fn|):Left=left+5:Top=top+5:Right=right-5:Bottom=bottom-5:MidM Bottom+Top)/2;x_unit_lcngth=(Right-Lcft)/(N-f-2);y_unit_length=( Bottom-Mid )/max;int gdriver = DETECT, gmode:initgraph(&gdriver. &gmode."”)：rectangle(left.top.right.botto!n)；setcolor( YELLOW):setcolor(GREEN):int x;x=Lcfl+x_unil_lcnglh/2;n=0:while(n<length) line(x.Mid.x.Mid-lnay_unit_length);x+=x_unit_length;n-M-;setcolor( LIGHTGRAY):getch();closegraphO;附录A3傅立叶变换有关算法本部分内容可以在 hg:2O2.3&75.33/dsD/C/PFT-FFT.cDD 获得。1、傅立叶变换数值计算函数/* 输入参数：f :击要进疔傅立叶变换的数值序列N :输入序列的长度M :频诰分析的频点数result:傅立叶变换的结果(蚁数序列)调用示例：M=100;DSFT(f0.N.MF);*/void DSFT(double N.int M,complex *result)int k.n;double omega.delta_omega:delta_onwga=2*yLPI/M;omega=0.0;for(k=0;k<M:k+)result(k=complex(0.0.0.0):for(n=0:n<N:n+)resultk+=f(n*exp(complex(0,-1 )*omega#n);oinega+=deka_omcga:2、卷积计算函数/ 参数说明fl:序列一：fl_bcgin:序列一开始下标：fl_end:序列一结束下标：f2:序列二：f2_begin:序列二开始下标：f2_end:序列二结束下标f :存放结果序列：begin: (f放结果序列开始卜标：end :存放结果序列结束下标调用示例：JuanJi(fl .0.10,f2.0.3.f3.&begin.&end)：/void JuanJi(double *fl,int fl_begin.ini fl_end.double *f2.ini f2_begin,int f2_end.double *f.int f_begin.intint m.n.i.begin.end:*f_begin=(fi_begin<f2_begin)?fl_begin:f2_begin：*f_end=*f_bcgin+fl_end-fl_begin+t2_end-f2_begin;for( i=0.n=*Lbegin:n<=*f_cnd:n+. i+) begin=(n-f2_end>fl_begin)?(n-f2_end):fl_begin;end=(n-f2_begin<fl_end)?(n-f2_begin):fl_end:ti|=0:for( m=begin :in<=end: in+) fli+=fln】-fl_begin*f2(nm-f2_bcgin;3、FFT变换函数/ 参数说明：xxx :需要进行变换的数组： N :长度：FFtorlFFT :进行正变换或反变换的指示参数调用示例：如果您在阅读过程中发现疏湄和错谋诸您尽快和编石取得联系 cxh中国科宁技术人7电了丁程巧信息科学系:coP>Ilghi“（）（），for( i= 1: i <=256:i+=2) xxxiJ=1.0;xxx1.3.57.is real partxxxi+l=0.0;/xxx24.6& is image partFFT(xxxJ2&1);*/#definc SXVAP(x.y) (double temp=(x):(x)=(y);(y)=temp:void FFT(double xxx(Jnt Mini FFTorlFFT)/*FFT 1 IFFT -1*/ double temp9tempntempi;double wtemp.wr.wpr.wpi.wi.theta;int n.mmaxmq jstep.p;n=N«l; q=l;for(p= 1 :pvn:p+=2)if(q>P)SWAP(xxxq.xxxp): SWAP(xxx|q+l J.xxx|p+1);)m=n»l:while(m>=2&&q >m) q=m:m»= 1:q+=m:mmax=2:while(n>mmax)istep=2*minax;theta=2.0*M_PI/(FFTorIFFT*mmax): wtemp=sin(double)(0.5*theta); wpr=-2.0*uiemp*wtemp;wpi=sin(double)theta):for( wr= 1.0.wi=O.O jn= 1 un<minax;m+=2) fortp=m;p<=n:p+=istep)q=p+mmax: tempr=wr*xxxq-wixxxq+l ; tempi=wr*xxxq+l +wi#xxxq: xxxq=xxxp| 2mpr: xxx(q+l |=xxx(p+l J-tempi;xxxp+=tempr;xxx(p+l J+=tempi:wr=(wten)p=wr)*wpr-wi#wpi+wr; wi=wi*wpr+wtemp*wpi+wi;inmax=istep;附录A4滤波器有关算法本部分内容可以在 hg:/2O2.3&75.33/dsD/C/Filtcr.cpD 获得。1、多项式的乘方展开函数/* 参数说明：A : 一维数纽，为多项式的系数数组M :多项式的阶次N :乘方次数B : 一维数纽，输出参数，放展开结果MNB:输也展开以后多项式的阶次/int Pnpe(double *A.int M.int N.double *B.int &MNB)double C|20|;int i.NK.k.j:MNB=M*N:for(i=l:i<20:i+) (Ci=O.O;B|i=O.O;if(N=O)B| 1 |=1.0;retum 0: for(i=l;i<=M+l;i+) BiJ=Ai;if(N=l) return 0:NK=M+1:fort i=2;i<=N:i+) for (j=l;j<=M+l:j+)for(k=l ;k<=NK:k-M-)Ck+j-l=C|k+j-l+Aj*Bk|;NK=M+NK:for(k=l;k<=:NKJc+)(B(k=Ck; Ck=0.0;Ireturn 0;I2、多项式乘积展开函数/ 参数说明：A : 一维数纽，存放第一个多项式的系数M :第一个多项式的阶次B :-维数纽，存放第二个女项式的系数N :第二个多项式的阶次C : 一维数纽，存放结果多项式的系数MN :结果娄项式的阶次/int Ypmp(double A.int MJouble *B.int N.double *C,int &MN)如果您在阅读过程中发现疏漏和错谋.i右您尽快和编者収得联系ncMoik cxh中国科学技术人宁电r r.WJjfH息科宁 广媒体通伫9 (Copyright 2000)int i.j;MN二M+N:for (i=l:i<=20:i+)Ci=0.0;fbr(i= 1 ;i<=M+1 ;i+)for (j=l ;j<=N+l :j+)Ci+j-l=Ci+j-l+Ai*BUJ；return 0;I3、有理分式替换函数/* 参数说明：C : 一维数纽，存放H(s)分“多项式系数NO : H(s)分母多项式阶次F1 : 一维数纽，存放变换式的分子多项式系数F2 : 一维数纽.存放变换式的分母務项式系数NF :变换式的阶次A: 一维数纽.变换结果的分子多项式系数B: 一维数纽，变换结果的分母多项式系数Nob : H(z)的阶次/int Bsf(double C.int NO.double 吓1 .double *F2jnt NF.double MA.doublc *B.inl &Nob) double D2O|.E2OG2O：int I.jiiNNDNNENG;Pnpe(F2.NF.NO.A.Nob):for (I=l:I<=20;I+)BI=0.0;for (1=1 :I<=NO+1:1+)Pnpe(FlNEii.D.NND):Pnpc(F2.NENO-ii.E.NNE);Ypmp( D.NND.E.NNE.GNG )：for(j= 1:j<=NG+1 ;j+) Bj=BLj |+C I fGj ：return 0：4. Butterworth滤波器阶数计算函数int BN (double alpha_p. double alpha_s double oiniga_p double omiga_s) double tkt2.t3.t4:tl = pow( 10.0.0 l*alpha_p)1.0：12 = pow( 10.0.0.1 *alpha_s) 1.0:13 = log 10(omiga_p/omiga_s);t4 = Iogl0(sqrt(tl/t2)/t3;return (int)ceil(t4);5、Chebyshev滤波器阶数计算函数int CN (double alphas. double alpha_p. double omiga.s double omiga_p)double tl.t2.t3.t4;11 = pow( 10.0.0.1 *alpha_s)-l .0:t2 = pow( i 0.0.0.1 *alpha_p)-1.0:l3 = acosh(omiga_s/omiga_p);t4 = acosh(sqrt(t l/t2)/t3;return (int)ceil(t4);6、模拟滤波器传递函数产生程序/参数说明:AP wpNH:通帯衰减：AS :阻帘衰减:通帶临界频率：WS :阴带临界频率:模拟滤波器的阶数:一维数纽，存放H(s)的分母多项式系数BC:维父数数纽，心放H的分母多项式的根:H(s)的分子常数#define SSH(x) log(x+sqrt(x*x+1.0)#define CCH(x) log(x+sqrt(x*x-1.0) #definc SHl(x) (exp(x)-cxp(-x)/2.0 样 define CHl(x) (exp(x)+exp(-x)/2.0 void bcg( float AR float AS.lloat WP.float WS, int &N .float *H.complex * B.tloat & C) (float ESEYNWCQXY.P:int i.kjl;char butter;complex BH2O.B22O:float P【=3 14159265;cout«Mbutteruash/cheybeshev(b/c)?M«MnM;cin»butter:if( tolower( butter)=,b,)YN=.5*logl0(pow(10.,0*AP)-l)/(pow(100l*AS)J)"logl0(WP/WS)：N=int( fabs( YN)i 0.9999999);WC=WP/(pow(pow( 10.01 *AP)-1 .X0.1/rtoat(2*N);C=pow(WC,N);如果您在阅读过程中发现疏漏和错谋.请您尽快和编者取得联系nclwoik cxh中国科宁技术人宁电f f.WJjfu息科学系 广媒体通佇一 (Copyright 2000)for(i=l;i<=N:i+)1141;P=PI*(II+0.5)/N+0.5): B|iJ=WC*complex(cos(Phsin(P);if(tolower(butter)=c,)ES=sqrt(pow( 10.,AS/10.)-!.);E=sqrt(pow( 10., AP/10.)-1.);YN=(CCH(ES/E)/(CCH(WS/WP):N=int(YN=fabs(double)YN)+0.9999999);WC=WP;C=pow( WC,N/( E*pow(2.N-l);Q=L/E;for(i=l:i<=N:i+)X=SSH(Q)/N;Y=float(2*i-1 尸PI/floal(2N);Bi=Yomplcx(SHl(X)"in(Y)CHl(X 尸 cos(Y)WC:Ifor (i=l:i<=20:i+)B 1| i=complex(0.0.0.0)：B2 i |=complex(0.0,0.0);B l2=complex( 1.0.0)：if(N!=l)(for(i=2:i<=N:i+)for (k=l:k<=i-l:k 卄)B2k+l=Bl|k-Blk+l*Bil；for(k= 1:k<=i+1: k+)Blk|=B2(kl；B2k=complex(0.090.0);for(i=l;i<=N+l:i+)HiJ=reaI(BlIi)/C; for(i=l:i<=N;i+)cout « w M« Bi) « MnM;for(i= 1 ：i<=N+1 ；i+)cout« ' M« H(i/H(N+1 « ”fT;7、BESSEL 函数Bessel函数的级数近似计算方法一 double Bcssel(double x)int i=Lk=i:double s=l.t=x/2:while(t > 1.0c-5)S4.=t*t；ii：k*=i;t=t*(x/2)/k; return s;Bessel数的级数近似计算方法:double bcssclO(doublc x)int i;double d.y,d2.sum;y = x/2.0:d = 1.0:sum = 1.0;for(i=l:i<=25:i 卄)d = d y/i;d2 = d * d:sum = sum + d2:if(d2 < sum% 1.0e-8) break:rctum(sum):Bessel函数的计算方法二 double mBesselF (double x)double ax. bessel:double y;if (ax=fabs(x) < 3.75) y = x/375;y J y;如果您在阅读过程中发现疏漏和借谋胡您尽快和编者取得联系nclwoik cxh中国科学技术大学电子工程与信息科学系垄媒体通信实验皇(Copyright 2000)bessel = 1.0+y*(3.5156229+y*(3.0899424+y*( 1.2067492 +y*(0.2659732+y*(0.360768e-1 +y 町.45813e2); else y = 3.75/ax;bessel = (exp(ax)/sqrt(ax)*(0.39894228+y*(0328592el+y WO.225319e-2+y(-0.157565e-2+y*(0.91628 le-2+y*(-0.2057706e-1 +y#(O.2635537e-l +y*(-0.1647633c-1+yP392377e2):Ireturn (bessel);备注:在VC中提供了0,l，.Jn和_yO,_yl,.,_yn尊系列函数用于计算两类Bessel函数,使用 方怯请自行察看联机帮助8、Kaiser 窗根据BESSEL负数得到Kaiser謝void kaiser(double s|.int n.int ni.double b)for(int i=l;i<n:i+)s(i|=Bassel(b*sqrt(l-(2.O*n/(m-l)>-l)*(2.O*n/(m-l )-1 )/Bassel(b);如果您在阅读过程中发现疏漏和错谋.胡您尽快和编者取得联系ncWk cxh