fftfrequency.c

对FFT的基2频析型算法实现

#include <ansi_c.h>
#include <math.h>
#include <cvirte.h>		/* Needed if linking in external compiler; harmless otherwise */
#include <userint.h>
#include "20041000.h"
#define PI 3.1415927

static int panelHandle;
int z,k,i,j;
double *a_re,*a_im,*b_re,*b_im,*c_re,*c_im;
int in;
int num,num1,num2;
int np; 
double u_re,u_im,w_re,w_im,t_re,t_im,p_re,q_re;
int nd2,np=0,nm1,m,Me,l,lpk,lp1,N;
int L;
double  d; 
int h;

int main (int argc, char *argv[])
{
	if (InitCVIRTE (0, argv, 0) == 0)	/* Needed if linking in external compiler; harmless otherwise */
		return -1;	/* out of memory */
	if ((panelHandle = LoadPanel (0, "20041000.uir", PANEL)) < 0)
		return -1;
	DisplayPanel (panelHandle);
	RunUserInterface ();
	return 0;
}

int CVICALLBACK set (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
         int m;
        
	switch (event)
		{
		case EVENT_COMMIT:
			 GetCtrlVal(panelHandle,PANEL_NUM,&num);
				if(num<0)
                          {
                          num=4;
                         SetCtrlVal(panelHandle,PANEL_NUM,num);
		          }
                        d=num;
			h=0;
                         while(d>1)
			{
			 d=d/2;
			 h++;
			}
			if(d<1.0)
			{
			m=1;
			for(k=0;k<h;k++)
			 m=m*2;
			 num=m;
			}
            SetCtrlVal(panelHandle,PANEL_NUM,num);    //将输入点数变换成2的整数次幂
            //c_re=malloc(num*sizeof(double));
            //c_im=malloc(num*sizeof(double));
            GetCtrlVal(panelHandle,PANEL_IN,&in);  
               switch (in)
	 	  {
		    case 0:
		         GetCtrlVal(panelHandle,PANEL_NUM,&num); 
		         num1=num;  
		         a_re=malloc(num1*sizeof(double));
		         a_im=malloc(num1*sizeof(double));
		    	 for(i=0;i<num1;i++)
			     {
			       a_re[i]=0;a_im[i]=0;
			     }
		    case 1:
		         GetCtrlVal(panelHandle,PANEL_NUM,&num); 
		         num2=num;  
		         b_re=malloc(num2*sizeof(double));
		         b_im=malloc(num2*sizeof(double));
		    	 for(j=0;j<num2;j++)
			     {
			       b_re[j]=0;b_im[j]=0;
			    }
			break;
		  }
			i=0;j=0;m=0;
                      break;
		}
	return 0;
}		   

int CVICALLBACK inall (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	switch (event)
		{
		case EVENT_COMMIT:
		switch (in)
		{
		    case 0:
		    	   GetCtrlVal(panelHandle,PANEL_REAL,&a_re[i]);
		    	   GetCtrlVal(panelHandle,PANEL_IMAGE,&a_im[i]);
		    	  i++;
		    	   SetCtrlVal(panelHandle,PANEL_REALNUM,i);
		    	   break;
		    case 1:
		    	   GetCtrlVal(panelHandle,PANEL_REAL,&b_re[j]);
		    	   GetCtrlVal(panelHandle,PANEL_IMAGE,&b_im[j]);
		    	  j++;
		    	   SetCtrlVal(panelHandle,PANEL_REALNUM,j);	
			break;
		}
			break;
		}
	return 0;
}

void efft(double *c_re,double *c_im,int N)
{				 
              	  nd2=N/2;
              	  nm1=N-1;
              	  m=1;
                      d=N;
			np=0;
                        while(d>1)
			{
			 d=d/2;
			 np++;
			}
          do{
              	  u_re=1;
              	  u_im=0;
             	   z=np+1-m;
              	  Me=pow(2,z);
             	   k=Me/2;
             	   w_re=cos(PI/k);
             	   w_im=-sin(PI/k);
             	   j=0;
           do{
            	    l=j;
	     	    do{
	        	   lpk=l+k;
	       	       t_re=c_re[l]+c_re[lpk];
	       	       t_im=c_im[l]+c_im[lpk];
	       	       q_re=c_re[lpk];
	        	   c_re[lpk]=(c_re[l]*u_re)-(c_re[lpk]*u_re)-(c_im[l]*u_im)+(c_im[lpk]*u_im) ;
	        	   c_im[lpk]=(c_re[l]*u_im)-(q_re*u_im)+(u_re*c_im[l])-(u_re*c_im[lpk]) ;
	         	   c_re[l]=t_re;
	        	   c_im[l]=t_im;
	        	   l=l+Me;
         	      }while(l<=N-1);
         	       p_re=u_re;
	        	   u_re=(u_re*w_re)-(u_im*w_im);
	        	   u_im=(p_re*w_im)+(w_re*u_im);
	       	       j=j+1;
	     	  }while(j<=k-1);
	       	      m=m+1;
	      }while(m<=np-1);
	       	      l=0;
	     	    do{lp1=l+1;
	         	    t_re=c_re[l]+c_re[lp1];
	          	    t_im=c_im[l]+c_im[lp1];
	         	    c_re[lp1]=c_re[l]-c_re[lp1];
	         	    c_im[lp1]=c_im[l]-c_im[lp1];
	         	    c_re[l]=t_re;
	         	    c_im[l]=t_im;
	         	    l=l+2;
	     	      }while(l<=N-1);
	             j=0,l=0;
	     	   do{
	     	          if(l>=j)
	           	        k=nd2;
	          	    else
	         	     {     	  
	         	     t_re=c_re[j];
	           	     t_im=c_im[j];
	          	     c_re[j]=c_re[l];
	          	     c_im[j]=c_im[l];
	          	     c_re[l]=t_re;
	          	     c_im[l]=t_im;
	           	     k=nd2;
	         	     }
	       	      while(k<j+1)
	       	         { 
	       	         j=j-k;
	           	     k=k/2;
	           	     }
	          	     j=j+k;
	          	     l=l+1;
	        	   } while(l<=nm1-1);
}

void eifft(double *c_re,double *c_im,int N)
{				 
     			  for(i=0;i<N;i++)
                  {
                    c_re[i]=c_re[i]/N;
                    c_im[i]=c_im[i]/N;
                  }				 
                  nd2=N/2;
              	  nm1=N-1;
              	  m=1;
                      d=N;
			np=0;
                        while(d>1)
			{
			 d=d/2;
			 np++;
			}
          do{
              	  u_re=1;
              	  u_im=0;
             	   z=np+1-m;
              	  Me=pow(2,z);
             	   k=Me/2;
             	   w_re=cos(PI/k);
             	   w_im=sin(PI/k);
             	   j=0;
           do{
            	    l=j;
	     	    do{
	        	   lpk=l+k;
	       	       t_re=c_re[l]+c_re[lpk];
	       	       t_im=c_im[l]+c_im[lpk];
	       	       q_re=c_re[lpk];
	        	   c_re[lpk]=(c_re[l]*u_re)-(c_re[lpk]*u_re)-(c_im[l]*u_im)+(c_im[lpk]*u_im) ;
	        	   c_im[lpk]=(c_re[l]*u_im)-(q_re*u_im)+(u_re*c_im[l])-(u_re*c_im[lpk]) ;
	         	   c_re[l]=t_re;
	        	   c_im[l]=t_im;
	        	   l=l+Me;
         	      }while(l<=N-1);
         	       p_re=u_re;
	        	   u_re=(u_re*w_re)-(u_im*w_im);
	        	   u_im=(p_re*w_im)+(w_re*u_im);
	       	       j=j+1;
	     	  }while(j<=k-1);
	       	      m=m+1;
	      }while(m<=np-1);
	       	      
	       	      l=0;
	     	    do{lp1=l+1;
	         	    t_re=c_re[l]+c_re[lp1];
	          	    t_im=c_im[l]+c_im[lp1];
	         	    c_re[lp1]=c_re[l]-c_re[lp1];
	         	    c_im[lp1]=c_im[l]-c_im[lp1];
	         	    c_re[l]=t_re;
	         	    c_im[l]=t_im;
	         	    l=l+2;
	     	      }while(l<=N-1);
	         	       j=0,l=0;
	     	   do{
	     	          if(l>=j)
	           	        k=nd2;
	          	  else
	         	     {     	  
	         	     t_re=c_re[j];
	           	     t_im=c_im[j];
	          	     c_re[j]=c_re[l];
	          	     c_im[j]=c_im[l];
	          	     c_re[l]=t_re;
	          	     c_im[l]=t_im;
	           	     k=nd2;
	         	     }
	       while(k<j+1)
	       	         { 
	       	         j=j-k;
	           	     k=k/2;
	           	     }
	          	     j=j+k;
	          	     l=l+1;
	        	   } while(l<=nm1-1);
}

int CVICALLBACK hfft (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	switch (event)
		{
		case EVENT_COMMIT:
		    DeleteGraphPlot (panelHandle, PANEL_GRAPH, -1, VAL_NO_DRAW); 
             if(in==0)
             {
             efft(a_re,a_im,num1); 
		     PlotY (panelHandle, PANEL_GRAPH, a_re, num1, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_RED);
		     PlotY (panelHandle, PANEL_GRAPH, a_im, num1, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_YELLOW);
		     }
		     else
		     {
		     efft(b_re,b_im,num2); 
		     PlotY (panelHandle, PANEL_GRAPH, b_re, num2, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_RED);
		     PlotY (panelHandle, PANEL_GRAPH, b_im, num2, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_YELLOW);
			 }
			break;
		}
	return 0;
}

int CVICALLBACK hifft (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	switch (event)
		{
		case EVENT_COMMIT:
		   DeleteGraphPlot (panelHandle, PANEL_GRAPH, -1, VAL_NO_DRAW); 
             if(in==0)
             {
             eifft(a_re,a_im,num1); 
		     PlotY (panelHandle, PANEL_GRAPH, a_re, num1, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_RED);
		     PlotY (panelHandle, PANEL_GRAPH, a_im, num1, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_YELLOW);
		     }
		     else
		     {
		     eifft(b_re,b_im,num2); 
		     PlotY (panelHandle, PANEL_GRAPH, b_re, num2, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_RED);
		     PlotY (panelHandle, PANEL_GRAPH, b_im, num2, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_YELLOW);
		     }
			break;
		}
	return 0;
}

int CVICALLBACK pzxs (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
    int o,p,q;
	switch (event)
		{
		case EVENT_COMMIT:
			 GetCtrlVal(panelHandle,PANEL_PUXIAN1,&o);
			 GetCtrlVal(panelHandle,PANEL_PUXIAN2,&p); 
			 GetCtrlVal(panelHandle,PANEL_PUXIAN3,&q);
	 		 switch(in)
	 		 {
	 		  case 0:
	 		  SetCtrlVal(panelHandle,PANEL_PXOUT1,sqrt(pow(a_re[o],2)+pow(a_im[o],2)));
	 		  SetCtrlVal(panelHandle,PANEL_PXOUT2,sqrt(pow(a_re[p],2)+pow(a_im[p],2)));
	 		  SetCtrlVal(panelHandle,PANEL_PXOUT3,sqrt(pow(a_re[q],2)+pow(a_im[q],2)));
	 		  break;
	 		  case 1:
	 		  SetCtrlVal(panelHandle,PANEL_PXOUT1,sqrt(pow(b_re[o],2)+pow(b_im[o],2)));
	 		  SetCtrlVal(panelHandle,PANEL_PXOUT2,sqrt(pow(b_re[p],2)+pow(b_im[p],2)));
	 		  SetCtrlVal(panelHandle,PANEL_PXOUT3,sqrt(pow(b_re[q],2)+pow(b_im[q],2)));
	 		  break;
	 		 }
			break;
		}
	return 0;
}
int CVICALLBACK ppxs (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)

{   
    double *p,*q;
	switch (event)
		{
		case EVENT_COMMIT:
		    switch(in)
	 		 {
	 		  case 0:
	 		  GetCtrlVal(panelHandle,PANEL_NUM,&num); 
		      num1=num;
		      p=malloc(num1*sizeof(double));
	 		  for(i=0;i<num1;i++)
	 		  p[i]=sqrt(pow(a_re[i],2)+pow(a_im[i],2));
	 		  PlotY (panelHandle, PANEL_GRAPH, p, num1, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_RED); 
	 		  break;
	 		  case 1:
	 		  GetCtrlVal(panelHandle,PANEL_NUM,&num); 
		      num2=num;
		      q=malloc(num2*sizeof(double));
	 		  for(j=0;j<num2;j++)
	 		  q[j]=sqrt(pow(b_re[j],2)+pow(b_im[j],2));
	 		  PlotY (panelHandle, PANEL_GRAPH, q, num2, VAL_DOUBLE,VAL_THIN_LINE, VAL_EMPTY_SQUARE, VAL_SOLID, 1, VAL_YELLOW); 
	 		  break;
	 		 }
			break;
		}
	return 0;
}

int CVICALLBACK jj (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{ 	
        double *t_re,*t_im,*r_re,*r_im;
        double *z_re,*z_im;
            d=num1+num2-1; 
			h=0;
            while(d>1)
			{
			 d=d/2;
			 h++;
			}
			if(d<1.0)
			{
			m=1;
			for(k=0;k<h;k++)
			 m=m*2;
			 L=m;
			}
 	        t_re=malloc(L*sizeof(double));
            t_im=malloc(L*sizeof(double));
            r_re=malloc(L*sizeof(double));
            r_im=malloc(L*sizeof(double));
            z_re=malloc(L*sizeof(double));
            z_im=malloc(L*sizeof(double));  
            switch (event)
	          	{
		      case EVENT_COMMIT:
                       for(i=0;i<L;i++)
                       {
                        t_re[i]=0; t_im[i]=0;
                        r_re[i]=0; r_im[i]=0;
                       }
                       for(i=0;i<num1;i++)
                       {
                        t_re[i]=a_re[i];t_im[i]=a_im[i];
                       }
                       for(j=0;j<num2;j++)
                       {
                        r_re[j]=b_re[j];r_im[j]=b_im[j];
                       }
                       efft(t_re,t_im,L);
                       efft(r_re,r_im,L);
                       for(k=0;k<L;k++)
                       {
                        z_re[k]=t_re[k]*r_re[k]-t_im[k]*r_im[k]; 
                        z_im[k]=t_re[k]*r_im[k]+r_re[k]*t_im[k];
                       }
                       eifft(z_re,z_im,L);
                       DeleteGraphPlot(panelHandle,PANEL_GRAPH,-1,VAL_IMMEDIATE_DRAW);
                       PlotY(panelHandle,PANEL_GRAPH,z_re,L,VAL_DOUBLE,VAL_THIN_LINE,VAL_EMPTY_SQUARE,VAL_SOLID,1,VAL_RED);
                       PlotY(panelHandle,PANEL_GRAPH,z_im,L,VAL_DOUBLE,VAL_THIN_LINE,VAL_EMPTY_SQUARE,VAL_SOLID,1,VAL_YELLOW);
	      	      break;
		        }
	             return 0;
}

int CVICALLBACK quit (int panel, int control, int event,
		void *callbackData, int eventData1, int eventData2)
{
	switch (event)
		{
		case EVENT_COMMIT:
			QuitUserInterface (0);
			break;
		}
	return 0;
}