devbox
我家自动化
这个屌丝很懒，什么也没留下！
关注作者
热门标签
热门文章
当前位置:   article > 正文

广义互相关求信号时延 JAVA实现_java广义互相关算法

作者：我家自动化 | 2024-06-25 12:36:11

赞

踩

java广义互相关算法

最近在做一个声音测量距离的系统，平台为Android，要用广义互相关来求信号的时延。

里面用到了FFT，FFT是参考这位仁兄的：

http://blog.csdn.net/syrchina/article/details/6670517

如何通过FFT求信号时延看这里：

http://zlgc.seu.edu.cn/jpkc2/ipkc/signal/new/course/two/2-6-11.htm


废话不多说了：上代码，下面的是FFT的

  1. public class FFT {
  2. 	private int N_FFT=0;				//傅里叶变换的点数
  3. 	private int M_of_N_FFT=0;			//蝶形运算的级数，N = 2^M
  4. 	private int Npart2_of_N_FFT=0;		//创建正弦函数表时取PI的1/2
  5. 	private int Npart4_of_N_FFT=0;		//创建正弦函数表时取PI的1/4
  6. 	private float SIN_TABLE_of_N_FFT [] = null;
  7. 	
  8. 	private static final double PI =  3.14159265358979323846264338327950288419716939937510;
  9. 	
  10. 	public FFT(int FFT_N)
  11. 	{
  12. 		Init_FFT(FFT_N);
  13. 	}
  14. 	public Complex[] complexLization(float data[])
  15. 	{
  16. 		Complex[] w = new Complex[data.length];
  17. 		for(int i=0;i<data.length;i++)
  18. 		{
  19. 			w[i].real = data[i];
  20. 			w[i].imag = 0.0f;
  21. 		}
  22. 		return w;
  23. 	}
  24. 	public float[] magnitude(Complex[] data)
  25. 	{
  26. 		float[] r = new float[data.length];
  27. 		for(int i=0;i<data.length;i++)
  28. 		{
  29. 			r[i] = (float) Math.sqrt(data[i].real * data[i].real + data[i].imag * data[i].imag);
  30. 		}
  31. 		return r;
  32. 	}
  33. 	//初始化FFT程序
  34. 	//N_FFT是 FFT的点数，必须是2的次方
  35. 	private void Init_FFT(int N_of_FFT)
  36. 	{
  37. 	    int i=0;
  38. 	    int temp_N_FFT=1;
  39. 	    N_FFT = N_of_FFT;					//傅里叶变换的点数 ，必须是 2的次方
  40. 	    M_of_N_FFT = 0;					//蝶形运算的级数，N = 2^M
  41. 	    for (i=0; temp_N_FFT<N_FFT; i++)
  42. 	    {
  43. 	        temp_N_FFT = 2*temp_N_FFT;
  44. 	        M_of_N_FFT++;
  45. 	    }
  46. 	    //printf("\n%d\n",M_of_N_FFT);
  47. 	    Npart2_of_N_FFT = N_FFT/2;		//创建正弦函数表时取PI的1/2
  48. 	    Npart4_of_N_FFT = N_FFT/4;		//创建正弦函数表时取PI的1/4
  49.  
  50. 	    //data_of_N_FFT = (ptr_complex_of_N_FFT)malloc(N_FFT * sizeof(complex_of_N_FFT));
  51. 	    //data_of_N_FFT = 
  52. 	    //ptr_complex_of_N_FFT SIN_TABLE_of_N_FFT=NULL;
  53. 	    //SIN_TABLE_of_N_FFT = (ElemType *)malloc((Npart4_of_N_FFT+1) * sizeof(ElemType));
  54. 	    SIN_TABLE_of_N_FFT = new float[Npart4_of_N_FFT+1];
  55. 	    CREATE_SIN_TABLE();				//创建正弦函数表
  56. 	}
  57. 	//创建正弦函数表
  58. 	private void CREATE_SIN_TABLE()
  59. 	{
  60. 	    int i=0;
  61. 	    for (i=0; i<=Npart4_of_N_FFT; i++)
  62. 	    {
  63. 	        SIN_TABLE_of_N_FFT[i] = (float) Math.sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);
  64. 	    }
  65. 	}
  66. 	
  67. 	private float Sin_find(float x)
  68. 	{
  69. 	    int i = (int)(N_FFT*x);
  70. 	    i = i>>1;
  71. 	    if (i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！
  72. 	    {
  73. 	        //不会超过N/2
  74. 	        i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);
  75. 	    }
  76. 	    return SIN_TABLE_of_N_FFT[i];
  77. 	}
  78. 	
  79. 	private float Cos_find(float x)
  80. 	{
  81. 	    int i = (int)(N_FFT*x);
  82. 	    i = i>>1;
  83. 	    if (i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！
  84. 	    {
  85. 	        //不会超过N/2
  86. 	        //i = Npart4 - i;
  87. 	        return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i];
  88. 	    }
  89. 	    else//i>Npart4 && i<N/2
  90. 	    {
  91. 	        //i = i - Npart4;
  92. 	        return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT];
  93. 	    }
  94. 	}
  95. 	private void ChangeSeat(Complex DataInput[])
  96. 	{
  97. 	    int nextValue,nextM,i,k,j=0;
  98. 	    Complex temp;
  99.  
  100. 	    nextValue=N_FFT/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法
  101. 	    nextM=N_FFT-1;
  102. 	    for (i=0;i<nextM;i++)
  103. 	    {
  104. 	        if (i<j)					//如果i<j,即进行变址
  105. 	        {
  106. 	            temp=DataInput[j];
  107. 	            DataInput[j]=DataInput[i];
  108. 	            DataInput[i]=temp;
  109. 	        }
  110. 	        k=nextValue;                //求j的下一个倒位序
  111. 	        while (k<=j)				//如果k<=j,表示j的最高位为1
  112. 	        {
  113. 	            j=j-k;					//把最高位变成0
  114. 	            k=k/2;					//k/2，比较次高位，依次类推，逐个比较，直到某个位为0
  115. 	        }
  116. 	        j=j+k;						//把0改为1
  117. 	    }
  118. 	}
  119. 	//FFT运算函数
  120. 	public void FFT(Complex []data)
  121. 	{
  122. 	    int L=0,B=0,J=0,K=0;
  123. 	    int step=0, KB=0;
  124. 	    //ElemType P=0;
  125. 	    float angle;
  126. 	    Complex W = new Complex();
  127. 	    Complex Temp_XX = new Complex();
  128.  
  129. 	    ChangeSeat(data);//变址
  130. 	    //CREATE_SIN_TABLE();
  131. 	    for (L=1; L<=M_of_N_FFT; L++)
  132. 	    {
  133. 	        step = 1<<L;//2^L
  134. 	        B = step>>1;//B=2^(L-1)
  135. 	        for (J=0; J<B; J++)
  136. 	        {
  137. 	            //P = (1<<(M-L))*J;//P=2^(M-L) *J
  138. 	            angle = (float)J/B;			//这里还可以优化
  139. 	            W.imag =  -Sin_find(angle);		//用C++该函数课声明为inline
  140. 	            W.real =   Cos_find(angle);		//用C++该函数课声明为inline
  141. 	            //W.real =  cos(angle*PI);
  142. 	            //W.imag = -sin(angle*PI);
  143. 	            for (K=J; K<N_FFT; K=K+step)
  144. 	            {
  145. 	                KB = K + B;
  146. 	                //Temp_XX = XX_complex(data[KB],W);
  147. 	                //用下面两行直接计算复数乘法，省去函数调用开销
  148. 	                Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag;
  149. 	                Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real;
  150.  
  151. 	                data[KB].real = data[K].real - Temp_XX.real;
  152. 	                data[KB].imag = data[K].imag - Temp_XX.imag;
  153.  
  154. 	                data[K].real = data[K].real + Temp_XX.real;
  155. 	                data[K].imag = data[K].imag + Temp_XX.imag;
  156. 	            }
  157. 	        }
  158. 	    }
  159. 	}
  160.  
  161. 	//IFFT运算函数
  162. 	public void IFFT(Complex []data)
  163. 	{
  164. 	    int L=0,B=0,J=0,K=0;
  165. 	    int step=0, KB=0;
  166. 	    //ElemType P=0;
  167. 	    float angle=0.0f;
  168. 	    Complex W = new Complex();
  169. 	    Complex Temp_XX = new Complex();
  170.  
  171. 	    ChangeSeat(data);//变址
  172. 	    //CREATE_SIN_TABLE();
  173. 	    for (L=1; L<=M_of_N_FFT; L++)
  174. 	    {
  175. 	        step = 1<<L;//2^L
  176. 	        B = step>>1;//B=2^(L-1)
  177. 	        for (J=0; J<B; J++)
  178. 	        {
  179. 	            //P = (1<<(M-L))*J;//P=2^(M-L) *J
  180. 	            angle = (float)J/B;			//这里还可以优化
  181.  
  182. 	            W.imag =   Sin_find(angle);		//用C++该函数课声明为inline
  183. 	            W.real =   Cos_find(angle);		//用C++该函数课声明为inline
  184. 	            //W.real =  cos(angle*PI);
  185. 	            //W.imag = -sin(angle*PI);
  186. 	            for (K=J; K<N_FFT; K=K+step)
  187. 	            {
  188. 	                KB = K + B;
  189. 	                //Temp_XX = XX_complex(data[KB],W);
  190. 	                //用下面两行直接计算复数乘法，省去函数调用开销
  191. 	                Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag;
  192. 	                Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real;
  193.  
  194. 	                data[KB].real = data[K].real - Temp_XX.real;
  195. 	                data[KB].imag = data[K].imag - Temp_XX.imag;
  196.  
  197. 	                data[K].real = data[K].real + Temp_XX.real;
  198. 	                data[K].imag = data[K].imag + Temp_XX.imag;
  199. 	            }
  200. 	        }
  201. 	    }
  202. 	}
  203. }
下面的是一个辅助类： 

  1. public class Complex {
  2. 	public float real;
  3. 	public float imag;
  4. 	
  5. 	public Complex()
  6. 	{
  7. 		this.real = 0;
  8. 		this.imag = 0;
  9. 	}
  10. }

下面是自相关求解过程源码： 

  1. public class Correlation {
  2. 		
  3. 	private Complex[] s1 = null;
  4. 	private Complex[] s2 = null;
  5. 	private int lag = 0;
  6. 	public Correlation(float sig1[] , float sig2[])
  7. 	{
  8. 		int maxLen = sig1.length > sig2.length ? sig1.length : sig2.length;
  9. 		maxLen = (int) (Math.log(maxLen)/Math.log(2.0)) + 1;//求出FFT的幂次
  10. 		maxLen = (int) Math.pow(2, maxLen);
  11. 		
  12. 		s1 = new Complex[maxLen];
  13. 		s2 = new Complex[maxLen];
  14. 		for(int i=0;i<maxLen;i++)
  15. 		{
  16. 			//这一步已经完成了补零的工作了
  17. 			s1[i] = new Complex();
  18. 			s2[i] = new Complex();
  19. 		}
  20. 		
  21. 		for(int i=0;i<sig1.length;i++)
  22. 		{
  23. 			s1[i].real = sig1[i];
  24. 			//System.out.println(s1[i].real);
  25. 		}
  26. 		for(int i=0;i<sig2.length;i++)
  27. 		{
  28. 			s2[i].real = sig2[i];
  29. 			//System.out.println(s2[i].real);
  30. 		}
  31. 		
  32. 		//求出信号的FFT
  33. 		float[] rr = new float[maxLen];
  34. 		FFT fft = new FFT(maxLen);
  35. 		fft.FFT(s1);
  36. 		fft.FFT(s2);
  37. 		conj(s2);
  38. 		mul(s1,s2);
  39. 		fft.IFFT(s1);
  40. 		
  41. 		float max = s1[maxLen-1].real;
  42. 		for(int i=0;i>maxLen;i++)
  43. 		{
  44. 			if(s1[i].real > max)
  45. 			{
  46. 				lag = i;
  47. 				max = s1[i].real;
  48. 			}
  49. 			//System.out.println(s1[i].real);
  50. 		}
  51. 	}
  52. 	public int getLag()
  53. 	{
  54. 		return lag;
  55. 	}
  56. 	//求两个复数的乘法，结果返回到第一个输入
  57. 	public void mul(Complex[] s1,Complex[] s2)
  58. 	{
  59. 		float temp11=0,temp12=0;
  60. 		float temp21=0,temp22=0;
  61. 		for(int i=0;i<s1.length;i++)
  62. 		{
  63. 			temp11 = s1[i].real ; temp12 = s1[i].imag;
  64. 			temp21 = s2[i].real ; temp22 = s2[i].imag;
  65. 			s1[i].real = temp11 * temp21 - temp12 * temp22;
  66. 			s1[i].imag = temp11 * temp22 + temp21 * temp12;
  67. 			//s1[i].real = s1[i].real * s2[i].real - s1[i].imag * s2[i].imag;
  68. 			//s1[i].imag = s1[i].real * s2[i].imag + s1[i].imag * s2[i].real;
  69. 		}
  70. 	}
  71. 	//求信号的共轭
  72. 	public void conj(Complex s[])
  73. 	{
  74. 		for(int i=0;i<s.length;i++)
  75. 		{
  76. 			s[i].imag = 0.0f - s[i].imag;
  77. 		}
  78. 	}
  79. }
下面的是测试代码： 

  1. public class Test {
  2.  
  3. 	public static void main(String args[])
  4. 	{
  5. 		//fft的测试函数代码
  6. 		/*int f0 = 10;
  7. 		int fs = 100;
  8. 		int FFT_LENGHT = 512;
  9. 		float[] testdata = new float[FFT_LENGHT];
  10. 		Complex[] r = new Complex[FFT_LENGHT*2];
  11. 		for(int i=0;i<r.length;i++)
  12. 		{
  13. 			r[i] = new Complex();
  14. 		}
  15. 		FFT fft = new FFT(FFT_LENGHT*2);
  16. 		
  17. 		for(int i=0;i<FFT_LENGHT;i++)
  18. 		{
  19. 			r[i].real = (float) Math.sin(2*3.1415926*i*f0/fs);
  20. 			r[i].imag = 0.0f;
  21. 			//testdata[i] = (float) Math.sin(2*3.1415926*i*f0/fs);
  22. 		}
  23. 		//r = fft.complexLization(testdata);
  24. 		fft.FFT(r);
  25. 		testdata = fft.magnitude(r);
  26. 		for(int i=0;i<FFT_LENGHT;i++)
  27. 		{
  28. 			System.out.println(testdata[i]);
  29. 		}*/
  30. 		
  31. 		//测试互相关函数
  32. 		
  33. 		int f0=10;
  34. 		int fs=150;
  35. 		int K = 50;
  36. 		int lag = 50;
  37. 		int dataLength = 16;
  38. 		float[] s1 = new float[dataLength];
  39. 		float[] s2 = new float[dataLength*2];
  40. 		
  41. 		for(int i=0;i<dataLength;i++)
  42. 		{
  43. 			//s1[i] = (float) Math.sin(2*3.14*(i*2.0/K/fs+f0)*i/fs);
  44. 			s1[i] = (float) Math.sin(2*3.14*f0*i/fs);
  45. 			//System.out.println(s1[i]);
  46. 		}
  47. 		for(int i=0;i<dataLength*1;i++)
  48. 		{
  49. 			s2[i] = 0;
  50. 		}
  51. 		for(int i=dataLength*1;i<dataLength*2;i++)
  52. 		{
  53. 			s2[i] = s1[i-dataLength*1];
  54. 		}
  55. 		Correlation correlation = new Correlation(s1, s2);
  56. 		//System.out.println("结果："+correlation.getLag());
  57. 	}
  58. }

再来个对应结果的图吧：

从图从左往右数，最高的刚好是第16个点，再看看测试代码里面的两个信号就知道刚好延迟了一个信号周期，刚好16个点，结果太完美了。

先写到这里，后面有 成果了再上传其它的代码。

本文内容由网友自发贡献，转载请注明出处：https://www.wpsshop.cn/w/我家自动化/article/detail/756200
推荐阅读
相关标签

Copyright © 2003-2013 www.wpsshop.cn 版权所有，并保留所有权利。

  

闽ICP备14008679号