c5_5_3.c

dsp入门与实践一书的源代码

#include "stdlib.h"
#include "stdio.h"
#include "math.h"
#define Pi 3.1416
#define N 8000
#define M 750

#pragma DATA_SECTION(SN,"bufferA")			
float SN[N];						
#pragma DATA_SECTION(Xf,"bufferB")			
float Xf[N]; 						
#pragma DATA_SECTION(x,"bufferE")			
float x[N]; 						
#pragma DATA_SECTION(y,"bufferF")			
float y[N]; 						
#pragma DATA_SECTION(result,"bufferC")			
int result[M]; 						

void main()
{
    int i,j,counter,counter2,Nu,Nu2,m,max,totalN,total2; 	
    int zpoint[M],zspace[M];
    int difzs[M],dd[M],ddd[M],result2[M];
    float k, f, a, b;						
    float fo, wo,sample,w1,w2,h,avg,rf;    
    float t,t1,t11,t3,t4,t2,t22;
    float k1[4],k2[4];    
    k = 0.3;    f = 0.52;    a = 0.05;    b = 0.1;           	
    fo = 60683.5;                				
    wo = 60000;                					
    sample = 1000000;          					
    w1 = 2*Pi*fo/sample;					
    w2 = 2*Pi*wo/sample;
    h = 2*Pi*fo/sample;
	srand((unsigned)time(NULL));				
    								
    for(i=1;i<N+1;i++)
    {
        Xf[i-1] = sin(w1*i);
        SN[i-1] = a*sin(w2*i);					
        t = (float)rand()/(float)RAND_MAX;			
        t = 2.0*t - 1;						
        SN[i-1] = SN[i-1] + b*t;				
    }       
    x[0] = 0.4;
    y[0] = 0.27;
    
    for (i=1; i<N; i++)
    {    	
    	k1[0] = y[i-1];
       	k2[0] = x[i-1] - powf(x[i-1],3) + f*Xf[i-1] - k*y[i-1] + SN[i-1];
        		for (j=0;j<2;j++)
        				{
        					t1 = y[i-1] + k2[j]*h/2;
	        				t11 = x[i-1] + k1[j]*h/2;
            				t3 = (Xf[i-1] + Xf[i])/2;
            				t4 = (SN[i-1] + SN[i])/2;
            				k1[j+1] = t1;
        					k2[j+1] = t11 - powf(t11,3) + f*t3 - k*t1 + t4;
         				}           
   	    	t2 = y[i-1] + k2[2]*h;
         	t22 = x[i-1] + k1[2]*h;
         	k1[3] = t2;
         k2[3] = t22 - powf(t22,3) + f*Xf[i] - k*t2 + SN[i];
         x[i] = x[i-1] + (k1[0] + 2*k1[1] + 2*k1[2] + k1[3])*h/6;
         y[i] = y[i-1] + (k2[0] + 2*k2[1] + 2*k2[2] + k2[3])*h/6;
     }     
     m = 0;

     for (i=0;i<N-1;i++)
     {
     	if ((x[i] >= -0.00000000001) && (x[i] <= 0.00000000001))
     	{
     	   	zpoint[m] = i;
        		m = m + 1;
        	}
       	else if (x[i]*x[i+1] < 0)
    		{	
    		 	zpoint[m] = i;
             	m = m+1;
        	}
     }     
     m = m - 2;

     for (i=0;i<m;i++)
          	zspace[i] = zpoint[i+1] - zpoint[i];
     
     m = m-1;
     for (i=0;i<m;i++)
          	difzs[i] = zspace[i+1] - zspace[i]; 
     for (i=0;i<m;i++)
          	difzs[i] = abs(difzs[i]); 
	i = 0;
	totalN = 0;

	while (i < m-1)
  	{
  		counter=0;
  		while (difzs[i] <= 3) 
     	{
     		if ( i > m-2)
     			break;
            	else
                i = i + 1;
        }          
        while (difzs[i] > 3)
        {
        	if (i > m-2)
         		break;
     		else
      			i = i + 1;
        }          
  		Nu = i;
  		while (difzs[i] <= 3)
  		{
  			counter = counter + 1;
    		if (counter > 15)
    			{
    			result[totalN] = zpoint[Nu];
       		totalN = totalN + 1;
       		break;
       		}
       		
       		if (i > m-2)
            		break;
            	else
                	i = i + 1;
       	}
     }     
     i = m - 1;
     total2 = 0;
     
     while (i > 0)
     {
     	counter2 = 0;
     	while (difzs[i] <= 3)
     	{
     		if (i<1)
     			break;
     		else
     			i = i - 1;
     	}     	
     	while (difzs[i] > 3)
     	{
     		if (i < 1)
     			break;
     		else
     			i = i - 1;
     	}     	
     	Nu2 = i;
     	
     	while (difzs[i] <= 3)
     	{
     		counter2 = counter2 + 1;
     		if (counter2 > 15)
     		{
     			result2[total2] = zpoint[Nu2];
     			total2 = total2 + 1;
     			break;
     		}     		
     		if (i < 1)
     			break;
     		else
     			i = i - 1;
     	}
     }          
     max = totalN;
     if (total2 > max)
     	max = total2;     	
     if ((max == 1)||(max == 0)||(max == 2))
     	avg = 0;
     else
     {
     	if (total2 > totalN)
     	{
     		totalN = total2;
     		for (i=0;i<total2;i++)
     		{
     			result[i] = result2[i];
     		}
     	}     	
     	for (i=0;i<totalN-1;i++)
     	     		dd[i] = abs(result[i+1]-result[i]); 
     	for (i=0;i<totalN-2;i++)
     	     		ddd[i]=abs(dd[i+1]-dd[i]); 
     	for (i=0;i<totalN-2;i++)
     	
     		if (ddd[i] > 250)
     		{
     			for (j=0;j<totalN;j++)
     					result[j]=0;
     			break;
     		}
         } 
       	avg = abs(result[totalN-1]-result[0]);
         avg = avg/(totalN-1);          
         rf = 2*Pi*avg/sample;
         rf = 2*Pi/rf;
         rf = fo + rf;
	}