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

1、中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&#

2、171; 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()r

3、eturn 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

4、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*( Com

5、plex &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

6、/l;y=y/l;return Complex(x.y);附录A2 BC下的绘图木部分内容可以右 htp:/2O3/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安装h

7、l录bgiegavga.bgi"文件支持如果您在阅读过程M1发现疏涓和iJi误.诘您尽快和编者取得联系nclwork cxh中国科学技术人宁电f息科学 广媒体通伫®、 (Copyright 2000)2、通用绘图函数void Plot(int leftjnt lop,ini 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我

8、示数纟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)；setcol

9、or( 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 :频诰分析的频点数r

10、esult:傅立叶变换的结果(蚁数序列)调用示例：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():for(n=0:n<N:n+)resultk+=f(n*exp(complex(0,-1 )*omega#n);oinega+=deka_omcga:2、卷积计算函数/ 参数说明fl:序列一：fl_bc

11、gin:序列一开始下标：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 m.n.i.begin.end:*f_begin=(f

12、i_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变

13、换函数/ 参数说明：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:

14、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&g

15、t;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;

16、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)do

17、uble 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、多项式乘积展

18、开函数/ 参数说明：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=

19、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.

20、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. d

21、ouble alpha_s double oiniga_p double omiga_s) double tkt2.t3.t4:tl = pow( 10.0.0 l*alpha_p)1.0：12 = pow( *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

22、tl.t2.t3.t4;11 = pow( *alpha_s)-l .0:t2 = pow( i *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

23、(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«

24、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息科学系 广媒体通佇一 (

25、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

26、;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()：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;

27、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(doub