fftrad2.asm

汇编语言编写的雷达动目标显示源程序

f14=f0*f6,  f12=f8+f12,                 f8=dm(i0,m2),   pm(i10,m10)=f9;
f13=f8+f12, f10=f8-f12,   			    f0=dm(i1,m1),   f7=pm(i8,m8);/*dm:sin*/
			f14=f11-f14,   			    dm(i2,m0)=f10,  f9=pm(i11,m12);
  /*start on next butterfly in each group*/
f12=f0*f7,  f3=f9+f14,  f9=f9-f14,      f6=dm(i0,m0),   f1=pm(i9,m9);/*pm:cos*/
f8=f1*f6,                               dm(i2,m2)=f13,  pm(i10,m10)=f4;
f11=f1*f7,                                              pm(i10,m10)=f9;
f14=f0*f6,  f12=f8+f12,                 f8=dm(i0,m0),   f7=pm(i8,m8);
end_group:
f12=f0*f7,  f13=f8+f12, f10=f8-f12,     f6=dm(i0,m0),   pm(i10,m13)=f3;

					
	    r4=r15+r2,          i1=b1;       /*PREPARE R4 FOR #OF BFLIES CALC*/
	    r15=ashift r4 by -1;             /*# OF BFLIES/GRP IN NEXT STAGE*/
	    r4=-r15,            i9=b9;
	    m1=r4;                                /*update inc for sin & cos */
	    m9=r4;
	    r5=ashift r5 by 1,  f1=dm(i1,m1);     /*update # bttrfly in a grp*/
	    r3=-r5;                               /*   inc for bttrfly branch*/
	    r13=r3-1,           m0=r3;            /*   wrap to 1st grp       */
	    r10=ashift r3 by 1, f7=pm(i9,m9);     /*   inc between grps      */
end_stage:  r9=r10-1,           m2=r13;           /*   wrap to grp +1        */

/*_________ next to last stage__________*/
	m1=-2;                  /*modifier to sine table pntr    */
	m8=r10;                 /*incr between groups            */
	m9=-2;                  /*modifier to cosine table pntr  */
	m10=r3;                 /*incr between bttrfly branches  */
	m12=r9;                 /*wrap around to next grp + 1    */
	m13=r13;                /*wrap to bgn of 1st group       */

	i0=refft+N-1;
	i1=sine+(N/2)-2;        /*pntr to 1st sine coeff         */
	i2=refft+N-1;
	i8=imfft+N-1;
	i9=cosine+(N/2)-2;      /*pntr to 1st cosine coeff       */       
	i10=imfft+N-1;
	i11=imfft+N-1;
										f0=dm(i1,m1),   f7=pm(i8,m8);
f12=f0*f7,                              f6=dm(i0,m0),   f1=pm(i9,m9);
f8=f1*f6,                                               modify(i11,m10);
f11=f1*f7,                                              f7=pm(i8,m12);
f14=f0*f6, f12=f8+f12,                  f8=dm(i0,m0);
f12=f0*f7, f13=f8+f12,  f10=f8-f12,     f6=dm(i0,m0);                           

/*Do the N/4 butterflies in the two groups of this stage*/

lcntr=N/4,  do end_group2 until lce;    
    f8=f1*f6,                   f14=f11-f14,    dm(i2,m0)=f10,  f9=pm(i11,m8);
    f11=f1*f7,  f3=f9+f14,      f9=f9-f14,      dm(i2,m0)=f13,  f1=pm(i9,m9);
    f14=f0*f6,  f12=f8+f12,                     f8=dm(i0,m2),   pm(i10,m10)=f9;
	f13=f8+f12, f10=f8-f12,  				    f0=dm(i1,m1),   f7=pm(i8,m8);
    f12=f0*f7,                  f14=f11-f14,    f6=dm(i0,m0),   f9=pm(i11,m12);

    f8=f1*f6,   f3=f9+f14,      f9=f9-f14,      dm(i2,m0)=f10,  pm(i10,m10)=f3;
    f11=f1*f7,                                  dm(i2,m2)=f13,  pm(i10,m10)=f9;
    f14=f0*f6, f12=f8+f12,                      f8=dm(i0,m0),   f7=pm(i8,m12);
end_group2:
    f12=f0*f7, f13=f8+f12,      f10=f8-f12,     f6=dm(i0,m0),   pm(i10,m13)=f3;


/*   The last stage      */


	i1=sine;        /*pntr to 1st sine coeff         */	
	i9=cosine;      /*pntr to 1st cosine coeff       */
	m2=N/2;                
	m0=(N/2)+1;		/* modifier */
	m4=1;
	m14=1;
	i0=refft+N/2;		/*pntr to REAL X0,X1*/
	i2=refft;			/*pntr to REAL X0',X1'*/
	i8=imfft+(N/2);		/*pntr to REAL X0,X1*/
	i10=imfft;			/*pntr to IMAG Y0',Y1'*/
	i11=imfft;  		/*pntr to IMAG Y0*/
	b3=refft;
	b4=refft;
	
	b12=imfft;			/* base register for 2nd butterfly*/
	b13=imfft;
	b14=imfft;	
	i3=refft+(3*N/4);	/*2nd bttrfly pntr to REAL X0,X1*/
	i4=refft+(N/4);		/*2nd bttrfly pntr to REAL X0,X1*/
	i12=imfft+(3*N/4);	/*2nd bttrfly pntr to REAL X0,X1*/
	i13=imfft+(N/4);	/*2nd bttrfly pntr to IMAG Y0',Y1'*/
	i14=imfft+(N/4);	/*2nd bttrfly pntr to IMAG Y0*/
	l3=N;
	l4=N;
	l12=N;
	l13=N;
	l14=N;


	/*start first bttrfly*/

f6=dm(i0,m2)/*X1*/		,f7=pm(i8,m14);/*Y1*/
f0=dm(i1,m4)/* sin */	,f1=pm(i9,m14);/* cos */
f8=f1*f6,/*CX1*/  m10=m0;
f12=f0*f7,/*SY1*/ m13=m2;
f11=f1*f7/*CY1*/ ,f12=f8+f12 /*CX1+SY1*/ ,f10=dm(i0,m0)/*X0*/ ,f9=pm(i11,m14);/*Y0*/
f14=f0*f6/*SX1*/ ,f6=dm(i3,m2)/*X1*/ ,f7=pm(i12,m14); /*Y1*/	
f8=f0*f6, f15=f11-f14;
f12=f1*f7, f13=f10+f12 ,f4=f10-f12;	
	
/*do two bttrflys in one loop with single twiddle fetch*/
/*use complex conjugated twiddel pairs N/4 appart	   */
lcntr=N/4 , do last_stage until lce;
f11=f0*f7, 	f3=f9+f15,	f2=f9-f15,		f10=dm(i3,m0),	f9=pm(i14,m14);
f14=f1*f6, 	f12=f8-f12, 				f6=dm(i0,m2),	f7=pm(i8,m14);
/*single twiddel fetch*/				f0=dm(i1,m4), 	f1=pm(i9,m14);
f8=f1*f6, 	f15=f11+f14, 				dm(i2,m2)=f13,	pm(i10,m13)=f3;
f12=f0*f7,	f4=f10+f12 ,f13=f10-f12,	dm(i2,m0)=f4,  	pm(i10,m10)=f2; 
f11=f1*f7,	f2=f9+f15,	f3=f9-f15,		f10=dm(i0,m0),	f9=pm(i11,m14);
f14=f0*f6,	f12=f8+f12,					f6=dm(i3,m2), 	f7=pm(i12,m14);
f8=f0*f6, 	f15=f11-f14, 				dm(i4,m2)=f13,	pm(i13,m13)=f3; 
last_stage: f12=f1*f7, f13=f10+f12 ,f4=f10-f12,		dm(i4,m0)=f4,	pm(i13,m10)=f2;




rts;    /*finished*/

mul:
	
	b2=refft;//回波信号实部的FFT
	l2=N;
	m2=1;
	b3=redata;// 回波信号实部数据
	
	
	b10=imfft;//回波信号的虚部的FFT
	l10=N;
	m10=1;
	b11=imdata;//回波信号虚部数据
	
	b1=recoff;//频域系数实部
	l1=N;
	m1=1;
	
	b9=imcoff;//频域系数虚部
	l9=N;
	m9=1;
	
	lcntr=N,  do end_mul until lce;//做1024点复乘
	f0=dm(i2,m2),f1=pm(i10,m10);//x1,y1//回波信号的实部FFT，虚部FFT
	f4=dm(i1,m1),f5=pm(i9,m9);	//x2,y2//系数的实部FFT，虚部FFT
	f9=f0*f4;   //复数乘法：实部乘实部-虚部乘虚部作为新的实部，实部乘虚部+虚部乘实部做为虚部。
	f12=f1*f5;
	f10=f0*f5;
	f13=f1*f4,f8=f9-f12;//新的实部
	f7=f10+f13;	//新的虚部		//f7->imag;f8->real
end_mul:	
	dm(i3,m2)=f7,pm(i11,m10)=f8;//将复乘结果放到回波信号的数组中替换掉并且将其实虚部交换，
	//后用新的回波信号中的数据做ＦＦＴ，为下一步的ＩFFT做准备。

	

	
	
	rts;
	//=================================除以1024=====================================//
ifft_last:
	
  	b2=refft;
//  	modify(i2,N-1);
//	l2=N;
	m2=1;
	b3=reifft;
	m3=1;
	
	b10=imfft;
//	modify(i10,N-1);
//	l10=N;
	m10=1;
	b11=imifft;
	m11=1;
	
	f11=2.0;
	lcntr=N*2,  do end_ifft until lce;
	f0=dm(i2,m2);               //取信号FFT实部放入f0中
	f12=1024.0;
	call divide;
	f15=f0;
	f0=pm(i10,m10);
	f12=1024.0;
	call divide;
	
end_ifft: dm(i3,m3)=f0,pm(i11,m11)=f15;
	
	rts;
//=====================================================================================//	
do_modual:
	
	
	
	f11=3.0;
	f1=0.5;
	lcntr=r0,do loop1 until lce;
		f2=dm(i2,m2),f6=pm(i10,m10);
		f3=f2;
		f7=f6;
		f2=f2*f3;
		f6=f7*f6;
		f13=f2+f6; //f13存平方和 
	    f4=rsqrts f13;     
             f12=f4*f4;
             f12=f12*f13;
             f4=f1*f4,f12=f11-f12;
             f4=f4*f12;                  //第一次迭代，相对误差为10-5
             f12=f4*f4;
             f12=f12*f13;
             f4=f1*f4,f12=f11-f12;
             f4=f4*f12;                  //第二次迭代，相对误差为10-7
             f12=f4*f4;
             f12=f12*f13;
             f4=f1*f4,f12=f11-f12;
		f4=f4*f12;						//第三次迭代，相对误差为10-7
		f4=f4*f13;
loop1:   dm(i7,m7)=f4;

		rts;
		
mti:
		b2=reifft;
		modify(i2,83);
		m2=1;
		b3=b2;
		modify(i3,543);
		m3=1;
	
		b10=imifft;
		modify(i10,83);
		m10=1;
		b11=b10;
		modify(i11,543);
		m11=1;
	
 		b1=remti;
 		m1=1;
 		
 		b9=imti;
 		m9=1;
 		
 		lcntr=460,do mti_loop until lce;
 		f1=dm(i2,m2),f2=pm(i10,m10);//结果相减。
 		f3=dm(i3,m3),f4=pm(i11,m11);
 		f3=f3-f1;    //相减后的实部
 		f4=f4-f2;     //相减后的虚部
mti_loop:	dm(i1,m1)=f3,pm(i9,m9)=f4;
 		

 divide:   f0=recips f12 , f7=f0;//recips为取倒数
          f12=f0*f12;
          f7=f0*f7,f0=f11-f12;       //第一次迭代，相对误差为10-3
          f12=f0*f12;             
          f7=f0*f7,f0=f11-f12;       //第二次迭代，相对误差为10-5 
          f12=f0*f12;                
          f7=f0*f7,f0=f11-f12;       //第三次迭代，相对误差为10-7 
          f0=f0*f7;      			  //无迭代时，仅保留第一条指令和这条指令
		  rts;		
.ENDSEG;