dft24.c

用于24点采样的傅立叶分解用的系数组,很适用的哦。这是我的一款电机保护中所采用的算法

void DFT_24(int *DFTDataAddr, float *DFT_Result)
/*以下是用于24点采样的傅立叶分解用的系数组，在程序中将他们已经分解处理。
  ConvCos[]={1024,989,886,724,512,265,0,-266,-512,-725,-887,-990,-1024,-990,-887,-725,-513,-266,-1,265,511,724,886,989};
  ConvSin[]={0,265,511,724,886,989,1024,989,886,724,512,265,0,-266,-512,-725,-887,-990,-1024,-990,-887,-725,-513,-266};
  
  按照DFT变换，有以下计算：
  Rx =989*{f(1)-f(11)-f(13)+f(23)}+
      866*{f(2)-f(10)-f(14)+f(22)}+
      724*{f(3)-f(9)-f(15)+f(21)}+
      512*{f(4)-f(8)-f(16)+f(20)}+
      265*{f(5)-f(7)-f(17)+f(19)}+
      1024*{f(0)-f(12)};
  Img=265*{f(1)+f(11)-f(13)-f(23)}+
      512*{f(2)+f(10)-f(14)-f(22)}+
      724*{f(3)+f(9)-f(15)-f(21)}+
      886*{f(4)+f(8)-f(16)-f(20)}+
      989*{f(5)+f(7)-f(17)-f(19)}+
      1024*{f(6)-f(18)};

  设：
      Tmp00 = f(0)-f(12);
      Tmp01 = f(1)-f(13);  Tmp02 = f(11)-f(23);
      Tmp11 = f(2)-f(14);  Tmp12 = f(10)-f(22);
      Tmp21 = f(3)-f(15);  Tmp22 = f( 9)-f(21);
      Tmp31 = f(4)-f(16);  Tmp32 = f( 8)-f(20);
      Tmp41 = f(5)-f(17);  Tmp42 = f( 7)-f(19);
      Tmp51 = f(6)-f(18);
  则：
      Rx =989*(Tmp01-Tmp02)+
          866*(Tmp11-Tmp12)+
          724*(Tmp21-Tmp22)+
          512*(Tmp31-Tmp32)+
          265*(Tmp41-Tmp42)+
          1024*Tmp00;
      Img=265*(Tmp01+Tmp02)+
          512*(Tmp11+Tmp12)+
          724*(Tmp21+Tmp22)+
          866*(Tmp31+Tmp32)+
          989*(Tmp41+Tmp42)+
          1024*Tmp51;				*/ 
{
	data union {long dlong;int dint[2]; char dchr[4];} DFTTmpData;
	int data Tmp1=0,Tmp2=0;    //函数内部的暂时变量
	int data DFT_R=0,DFT_X=0;  //用于放置DFT的实部和虚部
	
	//==  Rx:989*{f(1)-f(11)-f(13)+f(23)};   Img:265*{f(1)+f(11)-f(13)-f(23)}  =================
	Tmp1=*(DFTDataAddr+1)-*(DFTDataAddr+13);
	Tmp2=*(DFTDataAddr+11)-*(DFTDataAddr+23);
	
	DFTTmpData.dlong=990*((long)(Tmp1-Tmp2))>>10;
	DFT_R=DFTTmpData.dint[1];
	
	DFTTmpData.dlong=266*((long)(Tmp1+Tmp2))>>10;
	DFT_X=DFTTmpData.dint[1];
	
	//==  Rx: 724*{f(3)-f(9)-f(15)+f(21)};  Img:724*{f(3)+f(9)-f(15)-f(21)}  ===================
	Tmp1=*(DFTDataAddr+ 3)-*(DFTDataAddr+15);
	Tmp2=*(DFTDataAddr+ 9)-*(DFTDataAddr+21);
	
	DFTTmpData.dlong=725*(long)(Tmp1-Tmp2)>>10;
	DFT_R=DFT_R + DFTTmpData.dint[1];
	
	DFTTmpData.dlong=725*(long)(Tmp1+Tmp2)>>10;
	DFT_X=DFT_X + DFTTmpData.dint[1];
	
	//==  Rx :265*{f(5)-f(7)-f(17)+f(19)};  Img:989*{f(5)+f(7)-f(17)-f(19)};  ====================
	Tmp1=*(DFTDataAddr+ 5)-*(DFTDataAddr+17);
	Tmp2=*(DFTDataAddr+ 7)-*(DFTDataAddr+19);
	
	DFTTmpData.dlong=266*(long)(Tmp1-Tmp2)>>10;
	DFT_R=DFT_R + DFTTmpData.dint[1];
	
	DFTTmpData.dlong=990*(long)(Tmp1+Tmp2)>>10;
	DFT_X=DFT_X + DFTTmpData.dint[1];
	
	//==  Rx :886*{f(2)-f(10)-f(14)+f(22)};  Img:512*{f(2)+f(10)-f(14)-f(22)};  =================
	Tmp1=*(DFTDataAddr+ 2)-*(DFTDataAddr+14);
	Tmp2=*(DFTDataAddr+10)-*(DFTDataAddr+22);
	
	DFTTmpData.dlong=886*(long)(Tmp1-Tmp2)>>10;
	DFT_R=DFT_R + DFTTmpData.dint[1];
	
	DFTTmpData.dlong=512*(long)(Tmp1+Tmp2)>>10;
	DFT_X=DFT_X + DFTTmpData.dint[1];   
	
	//==  Rx :512*{f(4)-f(8)-f(16)+f(20)};  Img:886*{f(4)+f(8)-f(16)-f(20)};  ==================
	Tmp1=*(DFTDataAddr+ 4)-*(DFTDataAddr+16);
	Tmp2=*(DFTDataAddr+ 8)-*(DFTDataAddr+20);
	
	DFTTmpData.dlong=512*(long)(Tmp1-Tmp2)>>10;
	DFT_R=DFT_R + DFTTmpData.dint[1];
	
	DFTTmpData.dlong=886*(long)(Tmp1+Tmp2)>>10;
	DFT_X=DFT_X + DFTTmpData.dint[1];   
	
	//==  Rx :1024*{f(0)-f(12)};    Img:1024*{f(6)-f(18)};  =====================================
	Tmp1=*DFTDataAddr-*(DFTDataAddr+12);
	Tmp2=*(DFTDataAddr+ 6)-*(DFTDataAddr+18);
	
	
	
	*DFT_Result     =(float)(DFT_R+Tmp1);       //输出DFT的实部
	*(DFT_Result+1) =-(float)(DFT_X+Tmp2);      //输出DFT的虚部,
	                                          //由于傅立叶变换的参数正负号放反了，计算完毕后把它到过来
	                                          //2005年1月6日
}